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Abstract 

Unmanned  aerial  vehicles  (UAVs)  play  an  essential  role  for  the  US  Armed  Forces  by  per¬ 
forming  missions  deemed  as  “dull,  dirty,  and  dangerous”  for  a  pilot.  As  the  capability  of 
UAVs  expand,  they  will  perform  a  broader  range  of  missions  such  as  air-to-air  combat. 

The  focus  of  this  thesis  is  forming  trajectories  for  the  closing  phase  of  an  air-to-air 
combat  scenario.  A  UAV  should  close  with  the  suspected  aircraft  in  a  manner  that  allows 
a  ground  operator  to  visually  identify  the  suspected  aircraft  while  avoiding  visual/eletronic 
detection  from  the  other  pilot. 

This  thesis  applies  and  compares  three  methods  for  producing  trajectories  which  enable 
a  visual  identification.  The  first  approach  is  formulated  as  a  mixed  integer  linear  program¬ 
ming  problem  which  can  be  solved  in  real  time.  However,  there  are  limitations  to  the 
accuracy  of  a  radar  detection  model  formed  with  only  linear  equations,  which  might  justify 
using  a  nonlinear  programming  formulation.  With  this  approach  the  interceptor’s  radar 
cross  section  and  range  between  the  suspected  aircraft  and  interceptor  can  be  incorporated 
into  the  problem  formulation.  The  main  limitation  of  this  method  is  that  the  optimization 
software  might  not  be  able  to  reach  online  an  optimal  or  even  feasible  solution.  The  third 
applied  method  is  trajectory  interpolation.  In  this  approach,  trajectories  with  specified 
boundary  values  and  dynamics  are  formed  offline;  online,  the  method  interpolates  between 
the  given  trajectories  to  obtain  similar  maneuvers  with  different  initial  conditions  and  end- 
states.  With  this  method,  because  the  number  of  calculations  required  to  produce  a  feasible 
trajectory  is  known,  the  amount  of  time  to  calculate  a  trajectory  can  be  estimated. 
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Chapter  1 


Introduction 


Unmanned  aerial  vehicles  (UAVs)  such  as  the  Predator  have  proven  their  value  by  per¬ 
forming  reconnaissance  missions  and  accurately  launching  hcllhre  missiles  on  ground 
targets.  It  seems  feasible  that  in  the  near  future  UAVs  will  perform  additional  mis¬ 
sions  associated  with  fighter  pilots.  Such  a  vehicle  would  be  able  to  fight  aggressively 
and  if  shot  down,  would  not  have  some  of  the  undesirable  consequences  of  a  piloted 
aircraft.  For  instance,  a  search  and  rescue  mission  is  not  necessary  for  a  UAV. 

An  autonomous  UAV  which  can  perform  all  of  the  required  missions  with  profi¬ 
ciency  commensurate  to  a  human  pilot  would  be  very  difficult  to  design.  There  are 
many  situations  throughout  even  the  latest  conflicts  where  a  human’s  cognition  and 
situational  awareness  avoided  grave  consequences.  Major  Robert  Nolan  describes  one 
such  incident  during  the  first  Gulf  War  where  Capt  Landis  Cook,  an  A-10  pilot,  was 
given  clearance  to  fire  on  a  convoy  of  tanks  [21].  However,  Capt  Cook  decided  to 
visually  identify  the  targets  despite  the  aircraft  controllers’  assertion  that  there  were 
no  friendly  forces  in  the  area.  After  approaching  the  suspected  targets,  Capt  Cook 
identified  the  tanks  as  British  Challengers.  It  seems  that  only  a  UAV  with  the  trait 
of  what  Major  Nolan  refers  to  as  “judgment”  would  be  able  to  properly  handle  this 
situation.  The  current  level  of  technology  has  not  reached  this  difficult  goal,  and  as 
a  result,  a  human  in  the  loop  will  be  necessary. 

Including  a  human  in  the  decision  making  process  for  LTAVs  may  involve  more 
than  choosing  whether  or  not  to  destroy  a  target.  But,  certain  decisions  seemingly 
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can  be  made  within  the  cockpit  or  a  ground  control  unit.  Lt  Col  David  Brown’s 
experience  with  the  USAF’s  drone  program  for  four  and  a  half  years  during  the 
1990’s  demonstrates  this  point  [5].  While  flying  the  QF-106  from  the  ground  against 
manned  aircraft  in  a  dogfight  scenario,  Brown  made  the  observation,  “The  ONLY 
combat  capability  the  QF-106s  lacked  that  day  was  offensive  firepower!” 

UAV’s  with  the  capability  of  performing  an  air-to-air  intercept  could  have  been 
useful  in  recent  conflicts.  During  Operation  Desert  Storm,  the  Coalition  Forces  were 
able  to  quickly  obtain  and  maintain  air  superiority.  Coalition  aircrews  achieved  a  33-1 
air-to-air  kill  ratio,  and  Iraqi  pilots  essentially  ceased  air  combat  operations  within 
the  first  two  weeks  of  the  conflict  [10].  However,  as  many  as  121  Iraqi  aircraft  escaped 
to  Iran  when  Iraqi  pilots  realized  their  inability  to  fly  against  the  better  trained  and 
equipped  Coalition  pilots  [10]. 

Stopping  these  aircraft  from  fleeing  the  country  proved  to  be  a  very  difficult  task 
[10].  Coalition  fighter  pilots  guarded  the  border  between  Iraq  and  Iran  with  missions 
reaching  as  a  long  as  seven  hours.  Additionally,  airborne  warning  and  control  system 
(AWACS)  aircraft  were  unable  to  detect  all  of  the  Iraqi  aircraft  because  of  the  low 
altitudes  at  which  the  Iraqi  pilots  fled.  Finally,  Iraqi  pilots  would  often  wait  until  there 
were  holes  in  the  Coalition’s  barrier  or  when  there  was  no  barrier  at  all.  Considering 
the  limited  amount  of  Coalition  aircraft  and  pilots  which  may  have  been  assigned 
to  more  important  missions  such  as  close  air  support,  UAVs  capable  of  performing 
air-to-air  intercepts  would  have  been  ideal  for  this  situation. 


1.1  Objective 

A  UAV  may  be  able  to  conduct  successful  air-to-air  intercepts  even  without  the 
capability  to  dogfight.  During  the  Gulf  War  of  1991,  over  40%  of  the  air-to-air  kills 
were  performed  in  the  beyond  visual  range  region  mainly  due  to  improved  radar  and 
missile  technology  [20].  However,  if  dictated  by  the  rules  of  engagement  or  possibly  by 
an  imperfect  electronic  identification  system,  a  visual  identification  (VID)  is  necessary 
before  firing  a  missile.  A  VID  should  be  performed  in  minimum  time  while  reducing 
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the  chances  of  being  detected  by  the  suspected  aircraft. 

This  thesis  addresses  some  options  for  improving  a  fighter  pilot’s  ability  to  perform 
a  VID  on  a  single  suspected  aircraft.  A  computer  could  cue  the  pilot  to  an  optimal 
tra  jectory  which  could  decrease  the  amount  of  time  required  to  complete  the  maneuver 
and  reduce  the  chances  of  being  detected  by  radar  by  managing  range  and  attitude. 
Additionally,  this  scenario  could  be  adapted  to  a  UAV  equipped  with  a  camera  in 
order  that  a  ground  based  operator  could  ultimately  identify  the  suspected  aircraft. 

Three  methods  for  calculating  trajectories  which  enable  a  VID  are  applied  in  this 
thesis.  The  first  method  is  mixed  integer  linear  programming  (MILP),  which  has 
been  shown  to  be  a  powerful  and  effective  method  for  solving  rendezvous  problems 
in  real  time  [26].  The  drawback  of  this  approach  for  low- visibility  trajectories  is 
that  it  is  difficult  to  accommodate  certain  nonlinear  criteria,  such  as  radar  cross 
section,  involved  in  producing  a  stealthy  trajectory.  As  a  result,  the  problem  is  also 
formulated  as  a  nonlinear  program  using  a  direct  trajectory  optimization  method 
sometimes  referred  to  as  direct  collocation  [16].  Although  this  allows  a  more  accurate 
radar  avoidance  constraint,  a  real  time  solution  is  not  guaranteed.  A  trajectory 
interpolation  approach,  which  is  the  third  applied  method,  may  be  one  solution  to 
this  problem  [12].  A  library  of  maneuvers  parameterized  by  boundary  values  and 
aircraft  dynamics  are  calculated  offline;  online,  the  automatic  trajectory  synthesis 
interpolates  to  obtain  a  desired  trajectory. 

1.2  Thesis  Outline 

The  thesis  is  organized  as  follows:  Chapter  2  begins  with  a  description  of  the  different 
phases  of  air-to-air  combat.  It  describes  in  detail  the  different  tactics  used  for  the  clos¬ 
ing  phase  which  is  used  to  visually  identify  the  suspected  aircraft.  Chapter  3  presents 
a  brief  summary  of  the  different  techniques  used  for  trajectory  optimization  within 
the  past  35  years.  Chapters  4-6  contain  the  different  optimization  methods  used  in 
producing  trajectories.  Finally,  Chapter  7  concludes  the  thesis  with  a  summary  and 
suggestions  for  future  research. 
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Chapter  2 


Air-to-Air  Combat  Phases 

2.1  Five  Phases  of  Aerial  Engagements 

The  technology  and  tactics  of  air-to-air  combat  have  changed  significantly  since  the 
first  recorded  dogfight.  However,  Mike  Spick,  author  of  Fighter  Pilot  Tactics ,  states 
that  there  are  five  phases  of  aerial  engagements  which  have  remained  throughout 
air-to-air  combat’s  history:  (1)  detection,  (2)  closing,  (3)  attack,  (4)  maneuver,  (5) 
disengagement  [33] . 

Detection  is  extremely  important  for  air-to-air  combat  because  it  provides  the  pilot 
with  the  opportunity  to  make  the  first  decisions.  For  the  oldest  cases  of  aerial  en¬ 
gagements,  pilots  only  could  detect  the  enemy  visually.  Inexperienced  pilots  suffered 
the  highest  casualty  rates  by  focusing  on  flying  the  aircraft  rather  than  maintaining 
their  situational  awareness.  As  a  result,  many  pilots  were  shot  down  without  even 
seeing  their  adversaries. 

Detection  remained  essential  in  more  recent  versions  of  air-to-air  combat.  F4 
pilot  Bill  Jenkins  commented  about  the  importance  of  spotting  the  enemy  during  the 
Vietnam  war,  “The  MIG  21  is  a  very  small  airplane.  From  head  on  it’s  difficult  to 
see  at  more  than  2  miles  although  a  formation  is  easier.  American  aircraft  were  much 
larger  and  could  be  seen  at  greater  distances,  not  least  clue  to  the  smoky  exhaust 
trail”  [33].  Although  the  US  possessed  missiles  which  could  function  beyond  visual 
range  (BVR),  two  incidences  of  fratricide  forced  the  pilots  to  return  to  positive  visual 
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identification  [33]. 

Throughout  Operation  Desert  Storm  superior  technology  allowed  Coalition  fighter 
pilots  to  enjoy  significant  advantages  in  detecting  Iraqi  aircraft.  The  Airborne  Warn¬ 
ing  and  Control  System  (AWACS),  which  is  capable  of  detecting  aircraft  hundreds 
of  miles  away,  soon  became  know  as  the  “third  wingman”  [19].  Additionally,  non- 
cooperative  target  identification  systems  installed  on  aircraft  allowed  fighter  pilots 
to  determine  independent  of  an  AWACS  whether  an  aircraft  is  friendly  [19].  These 
advancements  directly  contributed  to  the  large  number  of  BVR  engagements.  How¬ 
ever,  William  Lewis  an  F-15E  pilot  states,  “This  ideal  long-range  engagement  cannot 
always  occur,  though,  because  of  enemy  defensive  maneuvering  and  the  constraints 
imposed  by  rules  of  engagement  [20].” 

Closing  is  the  second  phase  which  would  be  implemented  when  BVR  combat  is  not 
feasible.  A  beam  intercept  is  one  possible  maneuver  used  to  close  with  the  suspected 
aircraft.  This  method  decreases  the  chances  of  being  visually /electronically  detected 
by  the  enemy  while  increasing  the  chances  of  making  a  visual  identification  (VID). 

The  Israeli  Air  Force  completely  routed  Syria’s  fighter  planes  during  the  Bekaa 
Valley  conflict  in  June  of  1982  by  using  a  sound  strategy  for  the  closing  phase: 

Although  the  IAF  later  maintained  that  it  took  no  shots  at  Syrian  fighters 
from  beyond  visual  range,  it  evidently  made  extensive  use  of  blind-side 
tactics  by  employing  the  E-2C  to  vector  F-15s  and  F-16s  into  beam  attacks 
against  Syrian  MIGs,  where  their  radar  warning  systems  were  reportedly 
least  effective  [19]. 

At  the  end  of  the  conflict,  the  Israelis  had  destroyed  85  Syrian  aircraft  without  any 
losses  to  its  own  fighter  force  [19]. 

After  the  fighter  pilot  has  detected  the  target  and  begun  the  closing  phase,  he 
must  consider  when  to  fire  his  munitions,  which  is  referred  to  as  the  attack  phase. 
Throughout  World  War  I  and  much  of  World  War  II,  the  desired  attack  position 
was  from  behind  the  enemy’s  aircraft.  Attacks  from  astern  gave  the  fighter  pilot  the 
best  possible  firing  opportunity  as  well  as  a  good  defensive  position.  However,  this 
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position  was  no  longer  required  with  the  advent  of  homing  missiles.  Spick  now  defines 
the  attack  phase  as  “the  moment  when  the  attacker  starts  to  align  his  weapons  [33].” 

The  fourth  phase  is  maneuver;  this  phase  is  also  known  colloquially  as  “dogfight” 
[33].  Spick  states,  “It  is  the  most  spectacular  phase  of  combat  and  for  this  reason  is 
often  considered  to  be  the  most  important  [33].”  History  has  shown  that  this  is  not 
the  case;  only  about  one  in  five  victims  of  air-to-air  combat  are  shot  down  during 
the  maneuver  phase  [33].  Nevertheless,  maneuvering  is  crucial  when  the  first  three 
phases  do  not  finish  the  intercept.  This  would  be  the  most  difficult  phase  for  a  UAV 
to  accomplish  due  to  the  high  level  of  situational  awareness  required  to  win  a  dogfight. 

The  final  phase  is  disengagement,  which  although  important  for  a  UAV  is  not 
nearly  as  essential  as  for  a  manned  aircraft.  A  UAV  can  pursue  a  target  until  it 
expends  its  fuel,  whereas  a  piloted  aircraft  must  return  to  base  after  reaching  his 
“bingo”  fuel  state.  However,  if  the  UAV  must  leave  the  fight,  Spick  describes  two 
methods  of  disengagement: 

If  the  attacker  is  at  close  range,  i.e.,  800  yds  or  less,  the  hardest  possible 
turn  into  the  direction  of  attack  would  cause  the  enemy  fighter  to  over¬ 
shoot  it.  At  longer,  i.e.,  missile  ranges,  disengagement  can  be  made  by  a 
series  of  turns  hard  enough  to  make  it  difficult  for  the  attacker  to  achieve 
a  good  missile  bring  position,  but  not  so  hard  that  an  undue  amount  of 
speed  is  lost  [33]. 

The  five  phases  of  air-to-air  combat  parameterize  an  air-to-air  intercept.  This 
thesis  focuses  on  the  closing  phase  of  aerial  combat.  The  goal  is  to  produce  optimal 
trajectories  which  effect  a  VID  in  minimum  time,  while  reducing  the  chances  of  the 
target  detecting  the  interceptor. 

2.2  Closing  Phase 

As  previously  mentioned,  missile  technology  is  such  that  an  aircraft  can  be  destroyed 
in  the  BVR  region.  Merges  with  suspected  aircraft  will  still  occur  if  the  electronic 
identification  is  not  working,  or  the  rules  of  engagement  dictate  a  visual  identification. 
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In  order  to  visually  identify  the  aircraft,  a  UAV  could  carry  a  camera  slewed  to 
the  radar  track  of  the  suspected  aircraft.  This  will  allow  a  human  to  make  the  final 
decision  on  whether  or  not  to  destroy  the  target.  Ideally,  this  will  result  in  earlier 
VID’s  because  the  ground  operator  can  focus  on  identifying  the  aircraft  rather  than 
the  additional  task  of  flying. 

Four  tactics  for  performing  a  VID  are  described  in  this  section:  (1)  proportional 
navigation,  (2)  forward-quarter  intercept,  (3)  stern  conversion,  (4)  beam  intercept. 


2.2.1  Proportional  Navigation 


Using  proportional  navigation  is  potentially  the  fastest  method  for  an  interceptor  to 
merge  with  the  target  aircraft.  This  technique  is  commonly  used  for  missile  guidance 
but  can  also  be  used  by  an  aircraft.  The  objective  of  proportional  navigation  is  to 
minimize  |A|,  where  A  is  the  line  of  sight  angle,  shown  in  Figure  2-1. 


TARGET 


Figure  2-1:  The  top  figure  shows  an  example  target-interceptor  trajectory  formed 
using  proportional  navigation.  The  bottom  figure  shows  values  of  line  of  sight  angle 
A  at  specific  points  along  the  trajectory.  A  is  measured  as  the  angle  between  the 
x-axis  and  the  line  connecting  the  interceptor  and  target.  At  points  A  and  B  the 
interceptor  turns  in  order  to  obtain  a  constant  value  of  A  which  occurs  at  points  C 
and  D.  If  the  two  aircraft  continue  along  their  respective  headings  at  point  D,  they 
will  collide. 


The  guidance  law  for  two-dimensional  proportional  navigation  is  described  in  [38] 
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as 


(2.1) 


nc  =  N'VC\ 

where  nc  is  the  acceleration  command  which  is  perpendicular  to  the  interceptor’s 
velocity  vector,  N'  is  a  user  specified  gain,  Vc  is  the  the  closing  velocity,  and  A  is 
the  rate  of  change  of  the  line  of  sight  angle.  Vc  and  A  are  derived  with  a  few  simple 
equations  in  [38]. 


2.2.2  Forward-Quarter  Intercept 

When  a  visual  identification  is  required,  the  forward-quarter  (FQ)  intercept  is  also 
a  viable  tactic.  The  goal  of  this  maneuver  is  for  the  interceptor  to  position  itself 
with  a  good  firing  opportunity  while  denying  one  to  the  target.  The  FQ  intercept 
achieves  two  specific  objectives,  ft  obtains  a  certain  target  aspect  angle  (TAA),  which 
is  defined  in  Figure  2-2.  It  also  achieves  a  desired  amount  of  distance  to  be  traveled 


INTERCEPTOR 


/ 

/ 

/x 

TAA 

TARGET 

Figure  2-2:  Target  aspect  angle  (TAA)  is  the  angle  between  the  target’s  velocity 
vector  and  the  line  of  sight  vector  between  the  interceptor  and  the  target. 


along  the  final  collision  heading  which  is  from  point  2  to  point  3  in  Figure  2-3. 

There  are  several  advantages  to  using  this  tactic.  One  is  that  some  missiles  work 
more  effectively  from  an  FQ  firing  position  than  from  head-on.  Additionally,  an 
FQ  position  can  increase  the  range  of  visual  identification  because  more  area  of  the 
target  is  exposed  from  the  side  than  straight  on  [32],  However,  this  maneuver  may 
take  longer  to  achieve  a  firing  solution  than  if  the  aircraft  were  to  use  proportional 
navigation. 
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Figure  2-3:  The  FQ  intercept  consists  of  two  turns.  The  purpose  of  the  first  turn  at 
point  1  is  to  increase  the  TAA.  The  interceptor  makes  a  turn  away  from  the  target 
and  then  continues  along  a  straight  line  until  the  desired  TAA  is  obtained  at  point 
2.  Then,  the  interceptor  makes  another  turn  to  maintain  the  TAA.  Both  aircraft  are 
on  collision  course  at  point  3.  [32] 


2.2.3  Stern  Conversion 


Another  option  for  the  pilot  to  merge  with  a  suspected  aircraft  is  to  perform  a  stern 
conversion  [32],  This  maneuver  positions  the  interceptor  directly  behind  the  target 
aircraft,  which  is  a  highly  advantageous  position. 


1 


TARGET 

REQUIRED  { 
DISPLACEMENT 


Figure  2-4:  The  steps  to  determining  this  trajectory  occur  in  reverse  order  to  the 
actual  performance  of  the  maneuver.  The  first  is  to  decide  the  final  desired  range 
between  the  target  and  the  interceptor  once  the  maneuver  is  completed  at  point  4. 
The  second  step  is  to  determine  the  turn  radius  between  points  3  and  4  to  place  the 
interceptor  behind  the  target  at  the  desired  range.  Typically,  a  turn  which  points  the 
interceptor  at  the  target  throughout  the  majority  of  the  turn  is  chosen.  The  third 
step  is  to  calculate  the  conversion  range  and  required  displacement  at  point  3,  which 
will  enable  a  successful  completion  of  the  maneuver.  The  final  step  is  to  fly  to  the 
conversion  point,  which  consists  of  two  turns  occurring  at  points  1  and  2.  [32] 


The  major  advantage  to  performing  this  maneuver  is  that  the  target  will  have  great 
difficulty  obtaining  a  firing  opportunity  once  the  trajectory  is  completed.  However, 
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there  are  several  disadvantages  which  make  this  approach  atypical  for  pilots  to  use. 
This  is  by  far  the  longest  duration  trajectory  of  the  four  methods  for  performing  a 
VID.  Also,  with  the  advent  of  all  aspect  missiles,  it  is  unnecessary  for  the  interceptor 
to  fly  behind  the  target  to  obtain  a  firing  opportunity.  The  stern  conversion  is  still 
practiced  under  different  circumstances  such  as  refuelling  missions. 

2.2.4  Beam  Intercept 

In  order  to  visually  identify  the  suspected  aircraft,  a  fighter  pilot  will  typically  perform 
a  maneuver  called  a  beam  intercept  [11,  19,  24],  An  example  of  this  maneuver  is  shown 
in  the  Figure  2-5,  where  the  interceptor,  approaching  from  the  right,  is  performing 
the  maneuver  on  the  target,  on  the  left.  Three  regiments  characterize  this  trajectory 
at  the  final  state:  (1)  interceptor  is  a  prescribed  distance  away  from  the  target,  (2) 
interceptor’s  velocity  vector  is  pointed  at  the  target,  (3)  the  target  aspect  angle  (TAA) 
is  between  90  and  110  degrees  [11,  24].  For  the  simulations  presented  in  this  paper, 
the  desired  TAA  is  set  to  90°. 


Figure  2-5:  The  beam  intercept  usually  consists  of  two  turns.  The  first  is  a  turn  away 
from  the  target  to  obtain  the  desired  separation  range.  The  second  is  a  turn  towards 
the  target  to  point  the  interceptor’s  velocity  vector  at  the  target. 


There  are  several  advantages  to  performing  a  beam  intercept  as  opposed  to  using 
proportional  navigation,  which  might  allow  a  faster  merge: 
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■  Decreases  the  chances  of  being  detected  by  the  target’s  radar 

•  Potentially  decreases  time  to  VID,  because  looking  at  the  target  from  the  side 
rather  than  straight  on  will  increase  the  target  aircraft’s  visual  surface  area 

•  Increases  difficulty  for  the  target  to  visually  identify  the  interceptor,  because  during 
the  final  turn  of  the  maneuver  the  interceptor  is  essentially  pointed  at  the  target 

•  Reduces  the  closing  velocity  which  will  allow  more  time  for  the  VID,  without 
this  excess  time  the  interceptor  might  have  to  reposition  itself  to  obtain  a  bring 
opportunity 

•  Ultimately  places  the  interceptor  outside  of  the  target’s  weapons  cone 

•  Places  the  interceptor  in  an  offensive  position,  in  the  event  that  the  target  is  an 
enemy 

These  advantages  motivate  the  decision  in  this  thesis  to  use  the  beam  intercept  as 
the  maneuver  for  performing  a  visual  identification  on  a  suspected  aircraft. 

Radar  Model 

For  an  aircraft  flying  over  a  ground  based  radar,  it  is  possible  to  reduce  significantly 
the  amount  of  time  during  which  the  aircraft  may  be  detected  [23].  The  same  should 
be  true  in  performing  a  beam  intercept  with  a  possibility  of  avoiding  radar  detection 
completely.  There  are  instances  throughout  history  where  pilots  have  complained  of 
having  a  radar  lock  on  an  aircraft  and  then  losing  it  [34] .  In  some  cases,  this  was  the 
result  of  the  aircraft  changing  its  attitude,  which  altered  the  presented  radar  cross 
section  [34], 

The  amount  of  signal  energy  that  an  aircraft’s  radar  receives  is  approximated  in 
[34]  as 


signal  energy  =  k  av9 ^ — —  (2.2) 

where  k  =  (4^2),  Pavg  is  the  transmitter’s  average  power  output,  G  the  antenna 
gain,  <j  the  radar  cross  section  of  the  interceptor,  Ae  the  effective  area  of  the  antenna,  t 
the  time  the  radar  is  trained  on  the  interceptor,  and  R  the  range.  Of  these  parameters, 


the  interceptor  can  directly  affect  a  and  R  since  all  of  the  other  parameters  are 
dependent  on  the  target’s  radar  system.  This  equation  is  an  approximation  because 
the  true  amount  of  signal  energy  depends  on  the  radar  system’s  efficiency  in  analyzing 
the  data  [34],  When  integrating  the  radar  signal,  other  factors  such  as  noise  must  be 
taken  into  account  in  order  to  reduce  the  occurrence  of  false-positive  detections. 

In  order  to  avoid  radar  detection,  the  interceptor  can  affect  three  variables.  First, 
the  maximum  half  cone  angle  for  aircraft  radar  typically  ranges  between  45  to  60 
degrees.  As  a  result  the  interceptor  can  fly  a  significant  portion  of  the  beam  intercept 
outside  of  the  radar  cone  where  radar  detection  is  impossible.  The  second  and  third 
variables  which  the  interceptor  can  influence  are  range  and  radar  cross  section.  Within 
the  cone,  the  signal  energy  received  by  radar  is  inversely  proportional  to  range  RA 
and  directly  proportional  to  radar  cross  section  a 


signal  energy  cc 


a 


(2.3) 


In  performing  a  beam  intercept,  the  interceptor  can  maintain  a  large  distance 
between  itself  and  the  target  until  it  has  flown  outside  of  the  target’s  radar  cone. 
Additionally,  the  interceptor  can  affect  its  radar  cross  section,  which  is  a  function  of 
the  interceptor’s  azimuth  Kr  and  elevation  6r  angles  relative  to  the  target  aircraft. 
These  angles  are  determined  with  the  equations  in  [23] 
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(2.6) 


where  Rbe  is  a  rotational  matrix  which  rotates  the  unit  relative  position  vector  between 
the  target  and  the  interceptor  xe  from  an  earth  frame  to  a  body  frame  xb.  xb  is  a 
three  dimensional  vector  where  xb)i  represents  the  ith  component  of  xb.  Rbe  consists 
of  two  rotational  matrices  [23] 
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(2.7) 
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Lwe  =  sin  (f>  sin  7  cos  'I'  +  cos  0  sin  'I'  —  sin  0  sin  7  sin  'I'  +  cos  0  cos  T  sin  0  cos  7 
cos  0  sin  7  cos  \P  —  sin  0  sin  T  —  cos  0  sin  7  sin  T  —  sin  q i  cos  \P  cos  q i  cos  7 
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cos  a  cos  / 3  cos  a  sin  (3  —  sin  a 


Lbw  =  —  sin  (3  cos  f3  0 

sin  a  cos  f3  sin  a  sin  f3  cos  a 


(2.9) 


The  first  rotation  Lwe  is  from  the  earth  to  the  wind  axis  system,  and  the  second 
rotation  L^w  is  from  the  wind  to  the  body  axis  system  [23] .  The  angles  in  the  rotation 
matrices  represent  the  interceptor’s  attitude  where  7  is  flight  path  angle,  "0  heading, 
0  bank  angle,  a  angle  of  attack,  and  [3  sideslip.  The  body  fixed  coordinate  system 
is  defined  with  positive  x  axis  out  the  interceptor’s  nose,  positive  y  axis  out  the 
interceptor’s  left  wing,  and  positive  z  axis  follows  the  right  hand  rule  [23]. 

For  the  direct  collocation  development  described  in  chapter  5,  a  generic  aircraft 
radar  cross  section  model  is  used  similar  to  the  one  found  in  [23],  where  a  bivariate 
cubic  spline  is  used  to  map  the  interceptor’s  azimuth  and  elevation  angles  to  er,  as 
shown  in  Figure  2-6. 
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Figure  2-6:  Generic  aircraft  radar  cross  section  model  [23] 
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The  presented  radar  cross  section  can  vary  greatly  even  with  small  aircraft  vibra¬ 
tions  [34]  so  that  the  smooth  representation  shown  in  Figure  2-6  may  be  inaccurate. 
A  radar  cross  section  model  for  an  actual  aircraft  must  include  these  variations  to 
ensure  that  avoiding  radar  detection  is  still  possible. 

In  this  research,  it  is  assumed  that  the  interceptor  has  perfect  knowledge  of  the 
state  of  the  target.  This  information  could  be  provided  by  the  aircraft’s  own  radar 
or  an  AWACS. 
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Chapter  3 


Aircraft  Trajectory  Optimization 


There  has  been  a  great  deal  of  research  done  on  aircraft  trajectory  optimization  over 
the  past  35  years.  Initially,  the  calculus  of  variations  was  used  to  produce  optimal 
solutions  in  analytic  form,  which  are  solvable  in  real  time.  However,  solutions  of  this 
form  exist  for  very  few  problems,  and  the  models  describing  the  aircraft  dynamics 
need  to  be  simplified,  to  admit  analytic,  closed  form  solutions. 

If  an  analytic  solution  is  not  found,  the  problem  may  be  solved  using  a  shooting 
method  approach.  An  indirect  shooting  method  also  uses  the  calculus  of  variations 
to  cast  the  problem,  whereas  the  direct  shooting  method  does  not.  For  both  of  these 
methods,  the  iteration  process  in  finding  an  optimal  solution  involves  integrating  the 
trajectory  forwards  or  backwards  in  time  from  the  known  initial  states  to  the  desired 
final  states,  which  serve  as  boundary  conditions  to  the  solver. 

A  promising  method  for  trajectory  optimization  is  mixed  integer  linear  program¬ 
ming  (MILP).  For  this  approach,  the  aircraft  dynamics  and  constraints  are  modelled 
in  linear  form,  and  a  globally  optimal  solution  can  be  found  for  many  scenarios  in 
real  time. 

In  order  to  avoid  integrating  the  trajectories  as  required  with  shooting  methods, 
the  direct  collocation  approach  can  be  used.  With  this  method  the  problem  is  dis¬ 
cretized,  and  the  dynamics  of  the  system  are  approximated  with  a  set  of  algebraic 
constraints.  The  problem  is  then  cast  and  solved  as  a  nonlinear  program. 

The  final  method  considered  is  trajectory  interpolation.  This  new  approach  in- 
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volves  interpolating  between  two  similar  trajectories  to  obtain  a  family  of  maneuvers. 


3.1  Analytic  Solutions 

Analytic  solutions  for  optimal  control  problems  are  typically  derived  using  the  calculus 
of  variations,  which  forms  the  Euler  Lagrange  equations  [6].  An  optimal  solution  must 
satisfy  these  equations,  which  occur  in  the  form  of  a  two  point  boundary  value  problem 
(TPBVP)  where  some  boundary  conditions  are  given  at  the  initial  time  and  others 
at  the  final  time.  In  very  few  instances  an  optimal  solution  can  be  found  in  closed 
form.  Solutions  of  this  type  for  aircraft  trajectory  problems  typically  involve  flight 
constrained  to  the  horizontal  plane.  Research  in  this  area  mainly  consists  of  flight 
from  point  A  to  point  B  in  minimum  time  or  flight  from  one  heading  to  another. 

In  1971,  Erzberger  and  Lee  published  a  paper  describing  methods  for  producing 
three  types  of  trajectories  in  the  horizontal  plane  for  a  constant  speed  aircraft:  (1) 
point  A  to  point  B  with  designated  final  heading,  (2)  point  A  to  a  line  with  final 
heading  along  the  line,  (3)  point  A  to  point  B  with  unfixed  final  heading  [14],  The 
solution  for  each  of  these  problems  occurs  in  bang-off-bang  form,  where  the  aircraft  is 
either  flying  along  a  straight  path  or  turning  at  maximum  bank  angle.  Rather  than 
solving  a  TPBVP,  various  solutions  are  produced  geometrically  where  the  optimal 
solution  is  the  one  with  the  shortest  path  length.  The  authors  found  that  an  optimal 
trajectory  consists  of  no  more  than  four  turns. 

Clements  produced  in  1990  an  analytic  solution  for  the  problem  of  an  aircraft’s 
flight  in  the  presence  of  winds  from  point  A  to  point  B  with  final  heading  unfixed 
[8].  The  aircraft  flies  at  a  constant  altitude  and  constant  indicated  airspeed.  Using 
geometry  consisting  of  turning  circles  and  the  solution  to  a  fixed-point  equation,  the 
algorithm  is  able  to  produce  solutions  in  real  time.  To  follow  this  paper,  Clements 
also  solved  a  similar  problem  in  1992,  including  a  limit  on  maximum  roll  rate  [9]. 

The  problem  of  producing  trajectories  for  a  varying  speed  aircraft  from  point  A  to 
point  B  with  a  fixed  final  heading  has  also  been  addressed.  In  1994,  Ben-Asher  solved 
this  scenario  using  a  multiple  shooting  method  [2],  The  author  notes  the  challenges  of 
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onboard  implementation  due  to  the  difficulty  in  producing  reasonable  values  for  the 
costates  of  the  TPBVP.  To  follow  this  work,  in  1997  Shapira  and  Ben- Asher  produced 
analytic  solutions  to  a  similar  problem,  limiting  the  solution  space  to  bang-off-bang 
trajectories  and  posing  the  speed  of  the  aircraft  as  a  linear  function  of  the  aircraft’s 
turn  rate  [31].  Although  this  method  may  not  produce  the  true  optimal  solution,  a 
solution  is  guaranteed  in  real  time. 

Analytic  solutions  are  also  available  for  the  problem  of  changing  the  aircraft’s 
heading  to  a  specified  value  and  finishing  the  turn  with  a  desired  speed  [15,  35].  In 
1998,  Grimm  and  Hans  showed  that  their  method  for  various  scenarios  differs  only 
slightly  from  the  optimal  solution  obtained  using  a  direct  multiple  shooting  approach 
[15].  However,  in  order  to  maintain  an  analytic  solution,  a  constraint  on  load  factor 
was  not  included. 


3.2  Shooting  Methods 

Shooting  methods  are  commonly  used  to  solve  TPBVP’s,  and  consist  of  three  cate¬ 
gories:  direct  shooting,  indirect  shooting,  and  multiple  shooting  [3]. 

For  direct  shooting,  the  method  involves  solving  for  a  subset  of  the  initial  boundary 
values,  final  boundary  values,  and  parameters  which  create  an  optimal  trajectory. 
The  parameters  are  a  set  of  coefficients  for  the  controls,  which  must  be  written  in 
explicit  or  implicit  form.  For  each  iteration,  the  method  integrates  the  differential 
equations  of  the  system  forward  or  backward,  and  then  updates  the  variables  based 
on  the  resulting  objective  function  value  and  the  errors  in  the  boundary  conditions. 
This  method  works  especially  well  when  the  number  of  unknown  variables  is  small. 
However,  a  major  shortcoming  is  that  integrating  small  errors  in  the  initial  values 
can  lead  to  enormous  errors  at  the  endstate  [3]. 

With  the  indirect  shooting  method,  one  must  define  the  necessary  conditions 
for  optimality  in  the  form  of  Euler  Lagrange  equations.  The  set  of  values  to  be 
solved  are  the  costate  variables  at  the  initial  time  and  the  time  to  complete  the 
trajectory.  Often,  we  assume  a  form  for  the  propagation  of  the  costate,  and  the 
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terminal  state  constraints  are  mapped  onto  the  costate.  The  costate  is  then  solved 
while  the  states  are  propagated  forward  to  find  the  optimal  solution.  Typically,  it 
is  hard  to  provide  a  good  initial  guess  for  the  costate  variables  due  to  their  lack  of 
intuitiveness.  Additionally,  a  small  error  in  these  values  can  lead  to  highly  erratic 
values  at  the  endstate  [3]. 

In  order  to  address  the  instability  found  in  the  two  previous  methods,  a  multiple 
shooting  approach  was  introduced.  The  problem  is  discretized  into  multiple  time 
steps,  so  that  small  initial  errors  only  translate  from  one  time  step  to  the  next. 
Additional  constraints  are  added  in  order  to  ensure  continuity  of  the  states  at  each 
of  the  nodes.  The  problem  is  then  solved  as  many  TPBVP’s  using  either  the  indirect 
or  direct  approach  [3]. 

With  all  of  the  shooting  methods,  a  major  difficulty  is  incorporating  inequality 
constraints  that  are  a  function  of  the  states  along  the  trajectory  [3].  The  portion 
where  an  inequality  constraint  is  active  is  referred  to  as  a  constrained  arc.  In  order 
to  solve  these  types  of  problems,  the  number  and  location  of  constrained  arcs  must 
be  defined  beforehand  [6]. 

Shooting  methods  have  been  applied  to  many  different  aircraft  trajectory  prob¬ 
lems.  Seywald,  Cliff,  and  Well  in  1991  used  the  multiple  shooting  method  for  finding 
optimal  trajectories  for  an  aircraft  flying  in  the  vertical  plane  where  the  objective 
is  to  maximize  range  [30].  The  problem  is  constrained  with  maximum  thrust,  load 
factor,  and  dynamic  pressure  limits  where  boundary  conditions  specify  the  initial  and 
final  specific  energy,  altitude,  and  flight  path  angle.  The  dynamic  pressure  limit  adds 
constrained  and  unconstrained  arcs  to  the  problem.  The  authors  are  able  to  produce 
optimal  trajectories  for  various  scenarios,  but  they  mention  that  the  problem  has  not 
been  solved  for  “long  flight  times”,  due  to  the  difficulty  in  finding  an  appropriate 
switching  structure  for  the  costate  variables. 

In  1993,  Bocvarov,  Lutze,  and  Cliff  used  shooting  methods  to  determine  the  ef¬ 
fect  of  adding  thrust  vectoring  to  a  high  fidelity  F/A-18  model  in  performing  certain 
maneuvers  [4],  The  authors  note  difficultly  in  finding  a  good  initial  guess  for  the 
costate  variables,  the  time  to  complete  the  maneuver,  and  the  location  of  the  con- 
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strained  arcs.  In  order  to  address  this  problem,  the  authors  create  optimal  trajectories 
for  simpler  maneuvers.  Then,  the  boundary  conditions  or  parameters,  such  as  the 
amount  of  thrust  vectoring  power,  are  varied  slightly  towards  the  desired  values,  us¬ 
ing  the  previous  solution  as  an  initial  guess  for  the  costate  variables.  This  homotopy 
method  allows  an  efficient  means  to  calculating  trajectories  for  a  more  complex  set 
of  maneuvers  [3]. 

3.3  Mixed  Integer  Linear  Programming 

MILP  is  an  optimization-based  approach  which  can  be  used  to  solve  trajectory  gen¬ 
eration  problems  in  real  time.  The  problem  is  discretized  with  N  nodes,  where  linear 
equations  describe  the  vehicle’s  state  at  each  node.  Constraints  on  the  vehicle’s  states 
and  controls  can  be  added  as  long  as  they  occur  in  linear  form,  and  adjoint  binary 
variables  are  used  in  order  to  impose  or  relax  relevant  constraints. 

Recent  research  has  shown  MILP  to  be  a  promising  method  for  trajectory  opti¬ 
mization.  The  first  paper  using  this  approach,  published  in  2001,  outlines  the  al¬ 
gorithm  for  moving  multiple  vehicles  from  an  initial  state  to  a  desired  state  while 
minimizing  time  or  fuel  burn  [29].  The  method  ensures  that  the  vehicles  avoid  col¬ 
lisions  with  each  other  as  well  as  stationary  or  moving  obstacles.  Both  a  fixed  final 
time  and  a  receding  horizon  control  approach,  which  reduces  the  computation  time 
for  creating  paths,  were  tested. 

In  2002,  Richards  and  How  applied  the  previous  research  to  a  specific  vehicle 
model,  which  is  an  aircraft  with  a  maximum  speed  and  turn  rate  [25] .  The  aircraft  is 
also  constrained  to  a  constant  altitude.  In  this  paper,  the  authors  present  solutions 
for  path  planning  problems  where  multiple  aircraft  are  required  to  visit  waypoints 
while  avoiding  obstacles. 

In  order  to  reduce  the  computation  time  required  for  producing  aircraft  trajecto¬ 
ries  with  MILP,  a  new  receding  horizon  control  approach  was  introduced  in  2002  [1], 
Trajectories  are  calculated  using  a  combination  of  a  modified  version  of  Dijkstra’s  al¬ 
gorithm  in  addition  to  MILP.  This  method  has  the  advantage  of  avoiding  unescapable 
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obstacles,  and  it  is  also  a  candidate  for  producing  solutions  in  real  time  with  only  a 
small  reduction  in  optimality. 

Recently  in  2003,  Richards,  Kuwata,  and  How  successfully  demonstrated  that 
MILP  can  be  used  for  real-time  path  planning  [26].  The  authors  present  two  experi¬ 
ments  which  involve  small  ground  vehicles  about  the  size  of  a  remote  control  car.  The 
first  uses  receding  horizon  control  to  guide  a  vehicle  around  an  obstacle  to  a  desired 
final  state.  The  second  is  a  rendezvous  problem,  where  one  vehicle  travels  along  a 
straight  path,  and  the  other  maneuvers  to  reach  the  desired  relative  position. 


3.4  Direct  Collocation 

With  the  direct  collocation  method  for  trajectory  optimization,  the  states  are  approx¬ 
imated  with  cubic  polynomials,  and  the  dynamics  are  discretized  with  N  nodes  using 
collocation  of  the  states  and  controls  [16].  The  problem  is  formulated  as  a  nonlinear 
program  where  the  dynamics  are  approximated  with  algebraic  constraints. 

The  seminal  paper  describing  this  approach  was  published  by  Hargraves  and  Paris 
in  1987  [16].  They  demonstrated  the  method  with  three  numerical  examples,  one  of 
which  is  the  fighter  aircraft  minimum  time  to  climb  problem.  Their  optimal  trajectory 
for  this  maneuver  closely  matched  the  results  previously  obtained  using  shooting 
methods  [7]. 

In  1999,  Ringertz  demonstrated  the  performance  of  a  direct  collocation  solution 
with  live  flight  testing  [27].  A  pilot  flying  a  Saab  J35  Draken  administered  the  optimal 
commands  to  accelerate  the  aircraft  to  a  specified  airspeed;  the  test  results  were  very 
similar  to  the  computed  solution. 

To  follow  Ringertz’s  previous  work,  in  2000  he  produced  trajectory  solutions  for  a 
minimum  fuel  turn  [28].  The  aircraft  dynamics  are  approximated  as  a  reduced-order 
model  of  the  full  six  degrees  of  freedom  model  which  assumes  an  inertial  reference 
frame,  a  constant  aircraft  mass,  and  a  rigid  body.  The  optimal  trajectory  for  this 
maneuver  is  an  out  of  horizontal  plane  maneuver  where  the  aircraft  varies  its  altitude 
in  order  to  minimize  the  amount  of  fuel  burned.  The  author  discusses  the  difficulty 
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in  obtaining  solutions  in  real  time.  He  states  that  the  method’s  performance  can 
be  improved  through  solving  the  problem  with  a  small  number  of  nodes  and  then 
using  this  solution  as  an  initial  guess  for  solving  subsequent  problems  with  additional 
nodes. 

Horie  and  Conway  also  used  the  direct  collocation  technique  in  2000  to  design 
trajectories  for  an  evader  aircraft  to  reposition  itself  behind  a  pursuit  aircraft  [18]. 
This  tactic  consists  of  a  vertical  maneuver  using  post-stall  flight  and  thrust  vectoring 
to  slow  the  evader’s  velocity  and  allow  the  pursuer  to  pass  below.  Using  a  reduced- 
order  aircraft  model  similar  to  an  F-16’s  dynamics,  the  authors  are  able  to  produce 
solutions  with  various  boundary  conditions. 

In  2003,  Norsell  showed  that  it  is  possible  to  significantly  reduce  the  amount  of 
time  that  an  aircraft  is  detected  while  flying  over  ground-based  radar.  Using  the 
same  dynamic  model  as  Ringertz  in  [28],  the  author  produces  optimal  trajectories 
with  the  goal  of  decreasing  detection  time  and  minimizing  fuel  burn.  A  user  selected 
gain  varies  the  weight  on  the  importance  of  these  two  objectives.  For  this  type  of 
trajectory,  a  real  time  solution  is  not  essential  due  to  the  fixed  position  of  the  ground 
radar;  a  trajectory  can  be  calculated  offline  prior  to  the  actual  flight. 


3.5  Trajectory  Interpolation 

In  2004,  Dever,  Mettler,  Feron,  Popovic,  and  McConley  introduced  a  new  approach 
for  creating  trajectories  [12],  The  method  uses  sets  of  similar  maneuvers,  which 
have  different  boundary  conditions,  and  interpolates  between  them  to  obtain  a  set  of 
additional  trajectories.  The  new  paths  maintain  a  similar  analytic  structure  between 
desired  boundary  conditions;  a  desired  trajectory  with  specified  initial  conditions  and 
endstates  can  be  obtained  using  this  process.  For  example,  a  desired  trajectory  for 
a  helicopter  is  to  transition  from  an  initial  forward  velocity  to  a  hover  state.  Two 
trajectories  with  different  initial  forward  velocities  could  be  produced  offline;  online 
the  method  interpolates  between  them  to  obtain  a  feasible  maneuver  for  the  given 
initial  forward  velocity.  This  method  has  been  successfully  tested  on  a  three  degree 
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of  freedom  rotorcraft.  Trajectory  interpolation  potentially  has  the  major  advantage 
of  being  able  to  create  complex  maneuvers  in  real-time. 

3.6  Method  Selection 

Three  of  the  five  approaches  listed  above  are  tested  in  this  thesis.  It  is  unlikely 
that  an  analytic  solution  exists  for  creating  beam  intercept  trajectories  in  which  the 
interceptor  avoids  radar  detection.  Also,  due  to  the  inherent  difficulties  of  including 
constrained  arcs,  a  shooting  method  approach  is  not  included. 

The  methods  of  MILP,  direct  collocation,  and  trajectory  interpolation  are  applied 
in  this  thesis.  MILP  has  the  advantage  that  the  amount  of  time  to  converge  to 
a  solution  can  be  estimated.  A  shortcoming  of  this  approach  is  the  difficulty  in 
formulating  a  complex  aircraft  model  and  putting  constraints,  such  as  radar  detection 
avoidance,  into  a  linear  form. 

The  direct  collocation  method  has  the  advantage  of  ease  of  use.  The  problem  is 
simply  formulated  as  a  nonlinear  program  and  then  solved  with  optimization  software. 
However,  a  major  disadvantage  is  the  existence  of  multiple  local  minima.  The  method 
can  converge  to  a  solution  which  may  be  an  extremely  inefficient  trajectory.  These 
convergence  problems  make  an  onboard  implementation  difficult. 

The  trajectory  interpolation  approach  is  included  because  it  has  the  potential  to 
produce  trajectories  in  real  time.  Although  the  method  does  not  maintain  optimality, 
the  solutions  may  be  sufficiently  close.  For  onboard  implementation,  a  nearly  optimal 
beam  intercept  delivered  within  seconds  is  certainly  better  than  the  true  optimal 
trajectory  produced  a  day  later. 
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Chapter  4 


MILP  Approach 


In  this  chapter,  we  investigate  the  beam  intercept  problem  defined  in  Section  2.2.4 
using  mixed  integer  linear  programming  (MILP).  The  approach  closely  follows  the 
method  presented  in  [25]  with  variations  in  the  cost  function  and  constraints  to  reduce 
the  likelihood  of  the  target  detecting  the  interceptor.  While  the  rendezvous  problem 
of  [26]  models  the  interceptor  dynamics  in  a  reference  frame  relative  to  the  target,  the 
dynamics  for  this  approach  include  only  the  interceptor’s  states.  The  target  aircraft’s 
position  and  velocity  are  projected  to  produce  its  baseline  trajectory. 


4.1  Radar  Detection  Avoidance  Included  in  the 
Cost  Function 


Let  fXtk  and  fyj-  be  the  interceptor’s  applied  forces  in  the  inertial  x  and  y  directions 
at  node  k.  Solving  for  the  controls  fx  k  and  fy^  produces  the  trajectory  according  to 
the  following  dynamics  and  constraints.  Using  a  set  of  N  nodes,  the  dynamic  model 
of  the  aircraft  is  discretized  with  the  linear  equations 
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Figure  4-1:  Linear  constraints  are  used  to  approximate  the  maximum  speed  limit 
Vmax-  The  lines  m  —  1 ...  6  represent  each  of  the  linear  inequalities. 


where  (xint,k,  yint.,k)  are  the  interceptor’s  x  and  y  position  in  the  inertial  frame.  Simi¬ 
larly  (xint,k,yint,k)  and  ( Xint,k,yint,k )  are  the  interceptor’s  velocities  and  accelerations 
in  the  inertial  frame.  A t(k)  is  the  user  specified  time  between  each  node. 

The  maximum  allowable  speed  vmax  and  force  fmax  at  each  node  are  approximated 
with  the  linear  constraints 


7T 

<  Vmax  COS  —  (4.4) 

V  ke[l...N],me[l...M]  (4.5) 

A  fmax  COS  (4-6) 

V  k  e  [1...N-  l\,m  e  [1...M]  (4.7) 

where  M  is  the  number  of  linear  constraints  used  to  approximate  the  2-norm  on 
speed  and  thrust.  With  larger  values  of  M,  the  set  of  constraints  approaches  the 
actual  magnitude  constraints  on  velocity  and  force  vectors,  as  shown  in  Figure  4-1. 

In  the  MILP  approach,  the  operation  convenes  by  first  synthesizing  an  entire  set 
of  candidate  terminal  interceptor  states  at  each  point  in  time  termed  ( Xdes,k,ydes,k ), 
(xcies,k,  IJdes.k),  which  are  the  interceptor’s  position  and  velocity  relative  to  the  tar- 
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get’s  estimated  position  and  velocity  at  node  k.  Not  all  of  the  sets  ( Xdes,k,y<ies,k ), 
(xdes,k,ydes,k)  are  reachable  or  admissible.  The  search  process  involves  propagating 
the  interceptor  dynamics  such  that  they  eventually  approach  or  intersect  this  time 
history  of  terminal  points  {xdeStk,  ydeSjk),  ( xdeStk,ydeSjk )•  Any  intersection  thus  pro¬ 
duces  a  feasible  engagement.  The  least  expensive  of  all  these  feasible  trajectories 
that  intersect  the  set  of  candidate  terminal  states  is  output  as  the  trajectory. 

The  interceptor’s  desired  endstate  changes  in  order  to  include  a  maneuvering 
target,  whose  position  and  velocity  at  each  node  k  are  specified  before  the  MILP  step 


%des,k  %tgt,k  Rs  COS  (\& tgt,k  i  90  ) 

(4.8) 
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where  vs  is  the  final  desired  interceptor  speed.  (xtgt,  ytgt)  and  (xtgt,  ijtgt )  are  the  target’s 
x  and  y  position  and  x  and  y  velocity  respectively,  'i’tgt^k  is  the  target’s  heading  and 
Rs  is  the  desired  range  between  the  target  and  the  interceptor  at  the  end  of  the 
intercept.  The  interceptor’s  preferred  terminal  state  to  make  a  positive  identification 
is  ±90°  to  the  target  orientation,  which  is  shown  in  Figure  4-2.  Additional  binary 
variables  or  a  heuristic  can  be  used  to  select  which  desired  final  heading  and  position 
produces  the  best  solution. 
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Figure  4-2:  Potential  endstates  for  beam  intercept 


The  following  binary  constraints  model  the  completion  of  the  beam  intercept 
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where  Xint  represents  the  set  of  variables  {xint,k,  Vint,k),  {xint,k,  i/int,k)  and  Xdes  repre¬ 
sents  (. xdes,k,ydes,k ),  (ides^Vdes^)-  H  is  a  large  positive  number,  and  wk  is  a  binary 
variable  which  allows  the  constraints  to  be  relaxed.  This  formulation  admits  solutions 
that  drive  the  interceptor  to  intersect  the  candidate  set  of  desired  interceptor  states 
produced  earlier.  If  the  interceptor  has  reached  the  candidate  relative  position  and 
velocity  at  node  k,  then  wk  =  0  which  is  shown  in  Figure  4-3.  This  relaxation  of  the 
constraints  is  represented  in  the  objective  function  as  the  variable  tk.  At  node  k  if 
wk  =  0,  then  tk...N  =  0  otherwise  tk  =  1.  In  order  to  enforce  that  the  maneuver  is 
completed  at  some  node  k ,  an  additional  constraint  is  added 

N 

J2tk<N~  1 

fc=l 

The  optimization  objective  is  to  minimize  the  time  to  complete  the  beam  intercept 
and  the  accelerations  over  the  course  of  the  mission.  In  order  to  approximate  the 


(4.16) 
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Wj=0  w  =0  v\'n=0 

f  f 

[  Xinkk=l  Xinkk=2 

1  .  ^ 

Xtgt,k=l  Xtgt,k=2 

Figure  4-3:  This  figure  shows  the  candidate  endstates  of  the  interceptor  at  each  node 
k  for  a  target  with  a  constant  heading.  Xtntk  and  Xtgt)k  contain  the  states  of  the 
interceptor  and  target  respectively,  and  Rs  is  the  desired  separation  distance.  If 
there  is  a  feasible  solution  at  an  endstate,  the  binary  variable  wk  can  be  set  to  zero, 
otherwise  it  is  set  to  one  in  order  to  relax  the  desired  endstate  constraint.  MILP 
software  solves  multiple  linear  programming  problems  where  the  binary  variables  are 
varied  in  order  to  find  a  feasible,  optimal  solution. 


.f 

Xint,k=N 


Xtgt  ,k=N 


objective  of  maximizing  the  range  between  the  two  aircraft  within  the  linear  cost 
formulation  structure  that  MILP  permits,  the  term  c\yk\  is  introduced 

N 

min  J  =  y>|/x,fc|  +  e|/y|fc|  +tk-  c\yk\ )  (4.17) 

Jx,kify,k  . 

1=1 

Through  c\yk\,  the  interceptor  is  encouraged  to  maximize  its  cross- range  relative  to 
the  target  throughout  the  trajectory.  This  formulation  enables  a  penalty  for  proximity 
with  the  target.  The  target’s  initial  position  is  placed  at  the  origin  with  its  velocity 
vector  pointed  along  the  x-axis,  so  that  with  larger  values  of  c  the  interceptor  will 
maintain  higher  values  of  range  throughout  the  trajectory.  This  could  potentially 
decrease  the  chances  of  being  detected  by  the  target’s  radar  because  the  amount  of 
power  that  the  aircraft’s  radar  receives  is  heavily  dependent  on  range. 

In  Section  4.3,  the  problem  is  solved  once,  given  the  initial  states  of  the  interceptor 
and  the  target’s  projected  trajectory.  However,  the  problem  could  be  resolved  during 
the  maneuver  if  the  target  deviates  significantly  from  its  previously  calculated  path. 
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4.2  Additional  State  Constraint  for  Radar  Detec¬ 


tion  Avoidance 


Another  method  for  avoiding  radar  detection  involves  breaking  up  the  trajectory 
planning  problem  into  two  segments.  The  initial  portion  entails  that  the  interceptor 
maintain  at  least  a  certain  range  with  respect  to  the  target  until  it  crosses  the  target’s 
radar  cone.  We  mathematically  formulate  this  requirement  by  having  the  interceptor 
fly  to  a  point  at  least  a  certain  distance  RR  along  the  target’s  radar  cone,  as  shown  in 
Figure  4-4.  If  the  value  for  RR  is  sufficiently  large  given  the  radar  cross  section  of  the 
interceptor  and  the  capability  of  the  target’s  radar,  then  it  is  reasonable  to  assume 
that  the  interceptor  will  not  be  detected  by  the  target’s  radar.  Beyond  this  point,  the 
solution  for  the  intercept  problem  is  formulated  independent  of  range.  This  method 
has  the  advantage  over  the  previous  approach  in  that  the  interceptor  only  maintains 
distance  from  the  target  while  within  the  target’s  radar  cone.  Also,  it  is  easier  to 
chose  a  value  of  RR  than  the  coefficient  c  in  order  to  avoid  proximity  with  the  target. 

This  method  uses  the  equations  as  previously  described  but  with  an  additional 
requirement  that  the  interceptor  must  fly  to  two  waypoints.  The  c\yk\  term  is  also 
removed  from  the  cost  function.  The  first  point  is  a  distance  of  at  least  a  certain 
range  RR  along  the  target’s  radar  cone 


%tgt,k) 

—  RRcOS^RCA  +  \l/ tgt,k ) 

< 

Hvk 

(4.18) 

(%int,k  %tgt,k) 

+  RRcos(RCA  +  \l/ tgt,k ) 

< 

Hvk 

(4.19) 

( Vint,k  Vtgt,k  > 

1  T  RRsin^RCA  4-  ^tgt,k) 

< 

Hvk 

(4.20) 

V 

/c  e  [1 ...  IV] 

(4.21) 

where  RCA  is  the  radar  cone  angle  limit  and  vk  is  a  binary  variable.  Now,  tk...N  =  0 
if  both  v  and  w  are  zero  at  previous  time  steps.  The  second  point  is  the  previously 
defined  completion  state  of  the  maneuver. 
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TARGET 


INTERMEDIATE  WAYPOINT 


Figure  4-4:  In  order  to  reduce  the  chances  of  being  detected  by  the  target’s  radar, 
the  MILP  problem  is  formulated  with  an  additional  constraint  that  the  interceptor  is 
required  to  fly  to  an  intermediate  waypoint  of  at  least  a  distance  RR  along  the  radar 
cone.  RCA  indicates  the  maximum  angle  at  which  the  target’s  radar  will  function 
and  4/ tgt  is  the  target’s  heading. 


4.2.1  MILP  Software  and  Limitation  of  MILP  Approach 

Xpress  software  [36]  was  chosen  for  the  following  simulations  because  it  was  available 
at  Draper  Laboratory.  It  can  be  used  to  solve  linear,  quadratic  and  integer  program¬ 
ming  problems.  CPLEX,  another  MILP  optimization  software,  could  also  be  used  for 
this  approach. 

A  limit  on  minimum  speed  is  not  included  in  the  problem  formulation.  This  can 
result  in  infeasible  trajectories,  where  the  solution  produced  by  the  optimizer  may 
exceed  a  maximum  turn  rate  or  the  aircraft  might  stall.  A  minimum  speed  constraint 
with  the  same  accuracy  as  the  maximum  speed  approximation  would  require  N  x  Ad 
additional  binary  variables.  This  slowed  the  run  time  of  the  MILP  optimization 
software  Xpress  [36]  significantly.  Otherwise,  the  turn  rate  constraint  can  also  be 
satisfied  by  decreasing  fmax  [25]. 
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4.3  MILP  results 


For  the  following  the  simulations,  initially  the  target  and  interceptor  are  heading  to¬ 
wards  each  other  and  are  separated  by  27.5 km.  The  two  aircraft  have  the  same  initial 
speed,  and  the  target  is  turning  to  the  right  with  a  constant  turn  rate.  The  constants 
for  the  simulations  are  shown  in  Table  4.1.  Two  different  methods  using  MILP  are 
applied  to  solve  the  beam  intercept  problem.  The  first  involves  using  an  objective 
function  that  includes  time  to  completion  and  range  from  the  target  throughout  the 
maneuver.  The  second  MILP  approach  adds  an  additional  requirement  that  the  inter¬ 
ceptor  must  maintain  a  certain  range  from  the  target  until  it  crosses  the  radar  cone. 
The  objective  is  to  complete  this  task  and  finish  the  beam  intercept  in  minimum  time. 


Table  4.1:  Constants  for  MILP  beam  intercept  simulations 


N 

60 

M 

10 

mass 

10,442% 

H 

100,000 

V max 

167m/.s 

f max 

313  kN 

%int,k=l 

27.5  km 

yint,k=l 

0  km 

%tgt,k=l 

0  km 

Vtgt,k= 1 

0  km 

%int,k=l 

—  152.5  m/s 

yint,k= 1 

0  m/s 

%tgt,k=l 

152.5  m/s 

ijtgt,k= 1 

0  m/s 

Vs 

152.5  m/s 

Rs 

3.05  km 

e 

1  x  10"9 

A  t(k) 

2s 

4.3.1  MILP  Results  with  Adjusted  Cost  Function 

As  presented  in  Section  4.1,  the  cost  function  for  the  beam  intercept  scenario  is 

N 

min  J  =  y^(e|/x,fc|  +  e\fytk\  +  tk  -  c\yk\ )  (4.22) 

Jx,kify,k  .  A 

1=1 

where  the  c\yk\  term  encourages  the  interceptor  to  maintain  a  standoff  range  from 
the  target  vehicle. 

The  first  results  show  solutions  to  the  problem  with  variations  on  the  weighting 
factor  c  of  c\yk\.  Using  a  1.50  GHz  desktop  with  1.00  GB  of  RAM,  Xpress  software 
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(B) 


Figure  4-5:  Beam  intercept  using  MILP  with  y  weighting  coefficient  c  =  0.  Figure  (A) 
shows  the  interceptor  and  target  trajectory,  and  Figure  (B)  shows  the  interceptor’s 
controls  during  the  maneuver.  The  target  is  turning  at  a  constant  rate  of  0.3°/s.  A 
large  control  control  action  is  required  at  the  beginning  of  the  trajectory  in  order 
to  change  the  heading  of  the  interceptor.  For  the  majority  of  the  maneuver,  the 
interceptor  maintains  this  heading,  and  then  finally  a  large  control  effort  is  applied 
in  order  to  obtain  the  desired  relative  heading  and  position. 


[36]  is  typically  able  to  solve  this  problem  in  2-4  seconds. 

Figures  4-5  and  4-6  show  two  examples  of  a  beam  intercept  and  the  corresponding 
control  histories.  With  the  largest  value  of  the  y  weighting  coefficient,  the  maneuver 
requires  significantly  more  control  effort  and  the  time  to  completion  is  30  seconds 
longer  than  the  shortest  trajectory,  as  shown  in  Table  4.2.  In  both  scenarios,  at  the 
end  of  the  maneuver  the  interceptor  has  a  3.05 km  separation  distance  from  the  target 
and  a  terminal  heading  90°  different  than  target’s. 


49 


(A) 


time,  s 


Figure  4-6:  Beam  intercept  using  MILP  with  y  weighting  coefficient  c  =  0.1  which 
encourages  the  interceptor  to  increase  its  cross-range  from  the  target.  Figure  (A) 
shows  the  interceptor  and  target  trajectory,  and  Figure  (B)  shows  the  interceptor’s 
controls  during  the  maneuver.  The  target  is  turning  at  a  constant  rate  of  0.3°/s. 
Initially,  the  interceptor  makes  a  sharp  turn  to  the  right  in  order  to  obtain  a  larger 
distance  from  the  target.  Then,  the  interceptor  briefly  travels  along  a  straight  path 
and  banks  left  to  head  towards  the  target.  After  following  along  another  straight 
path,  a  small  turn  is  applied  to  complete  the  maneuver. 
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Figure  4-7:  With  higher  values  of  the  constant  c,  the  interceptor  maintains  a  larger 
distance  from  the  target  throughout  the  trajectory. 


Figure  4-7  shows  how  the  c\yk\  term  influences  the  range  between  the  interceptor 
and  target  during  the  beam  intercept.  The  four  curves  are  from  the  trajectories  in 
Table  4.2. 

A  range  of  trajectories  can  be  created  by  varying  c.  However,  there  is  no  guar¬ 
antee  that  the  interceptor  will  avoid  radar  detection.  After  the  trajectory  is  created, 
one  could  include  a  separate  evaluation  to  determine  if  the  interceptor  would  be  de¬ 
tected  by  the  target’s  radar.  If  the  trajectory  fails  the  test,  the  algorithm  could  run 
iteratively  with  larger  values  of  c. 


Table  4.2:  Results  from  MILP  beam  intercept  trajectories  including  c\y^\  term  in  cost 
function 


path 

length  (km) 

trajectory 
time  (s) 

computation 
time  (s) 

active  fx,fy 
duration  (s) 

min  /  max 
speed  (m/s) 

c  =  0 

14.2 

88 

2.9 

26 

150/159 

c  =  10  •  10“6 

14.7 

88 

4.7 

48 

153/167 

c  =  20  •  10"6 

15.0 

90 

5.2 

50 

153/167 

c  =  50  •  10“6 

17.8 

108 

5.3 

74 

150/167 

c  =  100  •  10“3 

18.2 

118 

3.4 

68 

106/167 
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4.3.2  MILP  Results  with  Additional  State  Constraint  for 


Radar  Detection  Avoidance 

For  this  formulation,  the  interceptor  is  required  to  fly  to  a  point  of  at  least  a  certain 
distance  RR  along  the  radar  cone,  where  the  radar  cone  angle  is  limited  to  45°.  The 
cost  function  for  this  scenario  is 


N 

min  J  =  YXe\fx,k\  +  e|/„)fcf  +  tk)  (4.23) 

Jx,kify,k  .  ^ 

1=1 

so  that  now  the  interceptor  must  fly  to  both  waypoints  in  minimum  time  with  a 
small  weight  on  control  effort.  With  the  additional  binary  variables  required  for  this 
approach,  the  computation  times  are  significantly  longer  as  seen  by  comparing  the 
results  in  Table  4.3  to  the  results  without  the  constraint  in  Table  4.2. 

Table  4.3:  Results  from  MILP  beam  intercept  trajectories  including  a  constraint  to 
fly  to  point  of  at  least  a  distance  RR  along  the  radar  cone  edge 


path 

length  (km) 

trajectory 
time  (s) 

computation 
time  (s) 

active  fx,fy 
duration  (s) 

min/max 
speed  (m/s) 

RR  =  9.2  km 

14.5 

90 

10.1 

32 

151/163 

RR  =  12.2  km 

16.2 

102 

29.9 

48 

124/165 

In  Figures  4-8  and  4-9,  the  interceptor  flies  to  a  point  of  exactly  the  required 
distance  along  the  radar  cone  and  finishes  the  trajectory  with  the  desired  relative 
heading  of  90°  and  range  of  3.05 km. 
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(A) 


(B) 


Figure  4-8:  Beam  intercept  using  MILP  where  interceptor  is  required  to  fly  to  a 
distance  of  at  least  9.2 km  along  the  target’s  radar  cone.  Figure  (A)  shows  the  air¬ 
crafts’  trajectories  along  with  the  radar  cone  limit  of  the  target,  and  Figure  (B)  is 
the  interceptor’s  control  history.  The  target  is  turning  at  a  constant  rate  of  0.3°/s. 
The  interceptor  begins  the  maneuver  with  a  right  turn  in  order  reach  the  end  of  the 
target’s  radar  cone.  Shortly  before  reaching  this  point,  a  hard  left  turn  is  executed  to 
wrap  around  the  radar  cone  and  to  point  the  interceptor  towards  the  desired  endstate. 
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Figure  4-9:  Beam  intercept  using  MILP  where  interceptor  is  required  to  fly  to  a 
distance  of  at  least  12.2 km  along  the  target’s  radar  cone.  Figure  (A)  shows  the 
aircrafts’  trajectories  along  with  the  radar  cone  limit  of  the  target,  and  Figure  (B)  is 
the  interceptor’s  control  history.  The  target  is  turning  at  a  constant  rate  of  0.3°/s. 
The  interceptor  performs  the  maneuver  in  a  similar  manner  as  Figure  4-8  but  with 
exaggerated  turns  in  order  to  fly  around  the  stronger  radar. 
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4.3.3  Discussion  of  Results 


As  stated  in  Section  2.2.4,  the  interceptor  can  avoid  radar  detection  from  the  target 
by  flying  outside  of  the  target’s  radar  cone,  maintaining  larger  distances  from  the 
target,  and  controlling  the  radar  cross  section  it  presents  to  the  target.  Both  MILP 
approaches  are  able  to  handle  the  first  two  criteria.  However,  they  are  limited  because 
it  is  difficult  to  include  the  radar  cross  section  of  the  interceptor  without  using  large 
amounts  of  binary  variables.  If  the  interceptor’s  maximum  radar  cross  section  is 
sufficiently  small,  this  value  could  serve  as  a  constant  in  setting  c  and  RR.  With 
larger  values  of  radar  cross  section,  a  method  which  directly  accounts  for  this  variable 
could  be  more  effective. 
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Chapter  5 


Direct  Collocation  Approach 

While  the  MILP  approach  can  create  feasible  trajectories  for  a  beam  intercept  in 
real  time,  it  is  difficult  to  formulate  a  linear  constraint  that  models  the  radar  cross 
section  of  an  aircraft.  The  following  direct  collocation  approach  can  address  this  issue, 
although  the  nonlinear  programming  formulation  might  result  in  longer  computation 
times  without  a  guarantee  of  converging  to  a  feasible  solution. 


5.1  Description  of  Direct  Collocation  Method 

With  direct  collocation  the  optimal  control  problem  is  cast  as  a  nonlinear  program 
(NLP)  with  constraints.  As  shown  in  Figure  5-1,  the  problem  is  discretized  with  N+ 1 
nodes  occurring  at  discrete  points  in  time,  (to,  where  t n  is  the  final  time 

tf- 

The  NLP  problem  is  written  in  terms  of  (N+l)*m  state  variables  consisting  of  x  = 
[xq,  x[,  x'2, . . . ,  x'N\  where  Xk  is  an  m-element  state  vector  defining  the  state  at  node 
k.  Additionally,  there  are  (IV +  1)  *  n  variables  for  the  controls  u  =  [u'0,  u\ ,  u2, . . . ,  u'N] 
where  Uk  is  an  n-element  column  vector  of  controls  at  node  k.  There  is  also  a  variable 
for  final  time  which  is  denoted  as  tf. 

The  objective  function  is  a  scalar  in  the  following  form 

J  =  u(x,  u,tf )  (5.1) 


57 


Figure  5-1:  Direct  collocation  involves  discretizing  the  trajectory  in  order  that  the 
states  and  controls  are  defined  at  specific  nodes.  The  nodes  are  not  constrained  to 
equal  spacing  in  time. 


where  u)  is  a  function  of  the  states,  controls,  and  final  time  and  is  subject  to  the 
following  constraints  of  system  dynamics,  boundary  conditions,  path  constraints,  and 
state  and  control  limits. 


The  dynamic  equations  describing  the  system  are  a  function  of  the  states  and 
controls 


X  =  f(x,u ) 

(5.2) 

There  are  also  constraints  describing  the  boundary  values  at  the 

initial  and  final  time 

bio  <  b(x0,u0,t0)  <bu0 

(5.3) 

bif  <  b(xN,uN,tN)  <  buf 

(5.4) 

Additionally,  the  solution  is  constrained  by  path  constraints 

9i  <  9(x,u,t)  <  gu 

(5.5) 

Finally,  there  are  limits  on  the  states  and  controls 

Xi  <  x  <  xu 

(5.6) 

Ul<U<Uu 

(5.7) 

The  boundary  conditions,  path  constraints,  and  limits  on  the  states  and  controls  are 


addressed  through  simply  entering  equations  (5.3-5. 7)  directly  as  they  appear  into 
the  NLP.  These  constraints  are  only  enforced  at  the  nodes,  and  the  trajectory  must 
be  checked  afterwards  to  ensure  feasibility.  If  the  trajectory  is  found  to  violate  a  path 
constraint,  additional  nodes  can  be  added  to  fix  the  problem. 

5.1.1  Collocation 

In  order  to  address  the  dynamics  of  the  system,  the  method  of  collocation  is  used 
to  approximate  the  equations  of  (5.2)  with  equality  constraints  [16].  The  first  step 
in  forming  these  constraints  is  to  represent  each  of  the  states  as  piecewise  cubic 
polynomials  between  the  nodes  in  the  form 

x  =  C0  +  C1S  +  C2S2  +  C3S3  (5.8) 

where  S  is  an  interpolation  variable  such  that  S  =  0  at  tk  and  S'  =  1  at  tk+i- 
Equation  5.8  and  its  derivative  are  evaluated  at  S  —  0,  S  —  1,  and  S  =  1/2  to  obtain 
the  following  equations 

xc  =  (xi  +x2)/2  +  T(f1  -  /2)/8  (5.9) 

x'c  =  -3(x1-x2)/2T-(/1  +  /2)/4  (5.10) 

A  =  fc  -  xc  (5.11) 

where  xc  is  the  interpolated  value  of  states  x\  and  x2  midway  between  the  nodes  at  t\ 
and  t2.  x'c  is  the  slope  of  equation  (5.8)  evaluated  at  S'  =  1/2.  /j,  f2,  and  fc  are  the 
same  as  equation  (5.2)  and  are  evaluated  using  (xi,Wi),  (. x2,u2 ),  (xc,uc)  respectively. 
uc  is  simply  Ul+U2 ,  and  the  variable  T  is  the  time  between  nodes. 

Equation  5.11  is  referred  to  as  the  defect,  and  a  set  of  N  *  m  of  these  equations 
enter  the  NLP.  The  states  and  controls  at  each  adjacent  node  are  varied  in  order 
to  set  the  defects  equal  to  zero.  Figure  5-2  shows  a  graphical  representation  of  this 
process  [16]. 

Once  the  NLP  program  has  been  solved,  a  continuous  control  history  can  be  found 
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Figure  5-2:  With  collocation  the  states  and  controls  at  adjacent  nodes,  for  example 
(xi,Ui)  and  (x2,u2),  are  varied  in  order  to  drive  the  defect  A  to  zero.  This  provides 
an  accurate  approximation  of  the  system  dynamics.  [16] 


by  linearly  interpolating  the  controls  at  adjacent  nodes  [16]. 

5.1.2  Scaling,  Analytic  Gradients,  and  Initial  Guess 

Optimization  software  typically  will  converge  more  efficiently  if  the  variables  are  of 
similar  magnitude.  For  the  simulations  in  this  chapter  the  states,  controls,  and  final 
time  are  scaled  with 


xs 

us 

tf. 


X 


Xmax 

U 


^ max 

tf 


(5.12) 

(5.13) 

(5.14) 


where  xmax,  umax ,  and  tmax  are  estimations  of  the  maximum  value  for  each  state, 
control,  and  final  time  respectively.  The  constraints  for  the  NLP  are  also  scaled 
considering  equations  (5.12-5.14)  in  order  to  maintain  consistency. 

In  the  following  simulations,  analytic  gradients  for  the  cost  function  and  all  con¬ 
straints  are  provided  to  the  optimization  software.  This  improved  the  convergence 
time  to  an  optimal  solution  significantly. 
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The  initial  guess  for  the  NLP  variables  is  set  to  half  of  the  estimated  maximum 
values  for  the  states,  controls,  and  final  time. 

5.2  Aircraft  Models 

Three  different  aircraft  dynamic  models  are  presented  for  this  approach:  (1)  the 
interceptor  maintains  a  constant  speed,  (2)  the  interceptor  varies  its  speed  while 
maintaining  a  constant  altitude,  (3)  the  interceptor  has  variable  speed  and  altitude. 
Simulations  for  the  simpler  models  typically  were  solved  with  less  computation  time; 
however,  the  extra  degrees  of  freedom  in  the  more  complex  models  allow  the  inter¬ 
ceptor  potentially  to  complete  the  maneuver  in  a  faster,  more  efficient  manner. 

5.2.1  Interceptor  Constant  Speed  Model 

The  dynamic  model  for  this  section  assumes  a  constant  speed  for  the  interceptor  and 
a  constant  speed  and  turn  rate  for  the  target.  Also,  the  interceptor  and  target  are 
at  the  same  altitude.  The  model  dynamics  are  described  in  a  relative  frame  with  the 
following  equations 


rel 

=  V  cos  T  -  Vtgt  COS  T tgt 

(5.15) 

Vrel 

=  v  sin  T  -  Vtgt  sin  T tgt 

(5.16) 

T 

tan(0)g 

V 

(5.17) 

*t9t 

=  d 

(5.18) 

where  xrei  and  yrei  are  the  x  and  y  position  of  the  interceptor  with  respect  to  the  tar¬ 
get,  v  and  vtgt  the  interceptor  and  target  speed  respectively,  0  the  interceptor’s  bank 
angle,  T  and  T tgt  the  interceptor  and  target  heading  respectively,  g  the  gravitational 
constant,  and  d  the  target’s  turn  rate.  The  origin  is  set  to  the  initial  starting  point  of 
the  target  with  \l/tgt(f0)  aligned  along  the  x-axis.  T  and  T tgt  are  defined  as  0°  along 
the  x-axis  and  increase  in  the  counter-clockwise  direction. 
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The  problem  is  formulated  as  an  NLP  with  the  following  variables: 

•  states  (m  =  3):  xrei(k),  yrei(k),  'i’(k) 

•  controls  (n  =  1):  (f>(k ) 

•  final  time:  tf 
where  Vfc  6  [0  . . .  N], 

The  objective  is  to  perform  the  maneuver  in  minimum  time  with  control  weighting 
to  assist  in  converging  to  a  locally  optimal,  feasible  solution 


N 

mm  J  =  if +pJ2(l)2(k)  (5-19) 

where  p  is  a  small  constant. 

The  constraints  for  the  NLP  include  dynamic  feasibility  which  is  approximated 
using  collocation  with  the  equations  in  (5.20).  These  defects,  shown  as  represent 
the  dynamics  of  equations  (5.15-5.17),  whereas  equation  (5.18)  is  easily  integrated  to 
find  the  target’s  heading  at  each  node.  An  example  of  calculating  these  equations  is 
shown  in  Appendix  A.  The  constraints  also  include  a  limit  on  the  maximum  bank 
angle,  radar  detection  avoidance,  and  final  desired  heading  and  position 


N) 


^k,j 

= 

0 

(5.20) 

V 

k  e  [l . . . 

N],Vj  G  [1  ...m\ 

4^min 

< 

VI 

S' 

4*  max 

(5.21) 

Rd{k ) 

< 

\JXrel(k ) 

+  Vrel(k ) 

(5.22) 

-V(N) 

= 

90° 

(5.23) 

X-rel  (  Af  ) 

= 

Rs  cos(T 

tgt(N)  +  90°) 

(5.24) 

Vrel(N) 

= 

Rs  sin('Li 

^(AO  +  90°) 

(5.25) 

where  Rd(k)  is  the  range  at  which  the  target’s  radar  can  detect  the  interceptor  at 
node  k,  and  Rs  is  the  final  separation  range  to  complete  the  beam  intercept.  Rd(k) 
is  defined  as 
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Figure  5-3:  (  is  the  distance  from  the  interceptor  to  the  edge  of  the  target’s  radar 
cone  and  is  defined  as  positive  (Figure  A)  if  the  interceptor  is  within  the  region  of 
detection  and  negative  otherwise  (Figure  B).  RCA  is  the  target’s  radar  cone  angle, 
and  vh tgt  is  the  target’s  heading. 


Rd(k) 
C  (k) 


ua 4 ( k ) 
l  _|_  e-<(k) 

cos(90°  -  ( RCA  +  ^tgt{k)))xrei(k)  - 
sin(90°  -  ( RCA  +  ^ tgt{k)))yrei{k) 


(5.26) 


(5.27) 


where  v  is  a  variable  dependent  on  the  power  of  the  target’s  radar,  cr(k )  the  inter¬ 
ceptor’s  radar  cross  section,  a  a  constant  which  is  explained  later,  ((k)  the  distance 
from  the  radar  cone  (shown  in  Figure  5-3),  and  RCA  the  radar  cone  angle. 

For  this  scenario,  a[k)  is  a  function  of  4>(k),  ^(/c),  xrei(k),  and  yrei(k)  where 
the  function  of  these  variables  is  described  in  Section  2.2.4.  It  is  assumed  that  the 
interceptor’s  angle  of  attack  is  nearly  a  constant  throughout  the  trajectory  and  does 
not  significantly  change  the  radar  cross  section  of  the  aircraft. 

With  the  detection  range  constraint,  the  sigmoid  form  of  equation  (5.26)  allows 
the  aircraft  to  fly  outside  of  the  radar  cone  without  concern  for  radar  cross  section 
or  range.  For  this  thesis,  a  =  16.39 /km  so  that  if  a  =  1  m2,  then  the  radar  detection 
range  can  be  described  with  the  curves  in  Figure  5-4  where  Rd  varies  according  to 
the  strength  of  the  target’s  radar.  With  larger  values  of  a,  the  radar  constraint 
becomes  more  like  a  binary  switch.  However,  if  this  value  was  set  too  large,  then  the 
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C  km 


Figure  5-4:  Radar  detection  range  as  a  function  of  distance  from  the  radar  cone  £ 
with  radar  cross  section  held  constant 

optimization  software  had  difficulty  in  converging  to  an  answer. 

5.2.2  Interceptor  Variable  Speed  Model 

The  time  to  perform  a  beam  intercept  can  be  reduced  significantly  if  the  interceptor 
is  permitted  to  vary  its  airspeed.  The  following  model  describes  the  dynamics  of 
the  system,  which  is  adjusted  from  [22]  in  order  to  formulate  the  problem  in  the 
relative  frame.  The  target  aircraft  is  assumed  to  have  a  constant  speed  and  turn  rate 
throughout  the  trajectory. 


%rel 

=  V  cos  T  -  Vtgt  COS  T tgt 

(5.28) 

Vrel 

=  v  sin  T  -  vtgt  sin 

(5.29) 

Tmax5  COS  (X  —  D 

V 

l-l 

(5.30) 

TmaxS  sin  a  sin  (p  +  L  sin  (p 

(5.31) 

/ IV 

^tgt 

=  d 

(5.32) 

64 


where  Tmax  is  the  interceptor’s  maximum  thrust  available,  5  the  throttle  setting,  a 
the  interceptor’s  angle  of  attack,  p  the  interceptor’s  mass,  and  L  and  D  are  the 
interceptor’s  lift  and  drag  respectively  and  are  calculated  with 

L  =  CLl-pv2S  (5.33) 

D  =  CDl-pv2S  (5.34) 

A  model  similar  to  an  F-16  found  in  [18]  is  used  for  the  wing  planform  area  S,  and 

the  lift  and  drag  coefficients 

CL  =  0.0174  + 4.3329a-  1.3048a2  +  2.2442a3 
— 5.8517a4(0  <  a  <  tt/6) 

CD  =  0.0476  -  0.1462a  +  0.0491a2  +  12.8046a3 
—  12.6985a4(0  <  a  <  tt/6) 

The  density  is  modelled  with  [18] 

p  =  ps[ 1  -  0.00688(/r/1000)]4'256  (5.37) 

where  h  is  altitude  measured  in  feet,  and  ps  is  the  density  at  sea  level. 

The  maximum  thrust  available  Tmax  is  a  constant  equal  to  the  weight  of  the 
interceptor  [18].  It  is  assumed  that  with  changes  in  airspeed  and  density,  the  true 
thrust  available  is  always  greater  than  or  equal  to  Tmax. 

Now,  the  nonlinear  program  has  the  following  variables: 

•  states  (m  =  4):  xrei(k),  yrei(k),  v(k),  'k(fc) 

•  controls  (n  =  3):  (f>(k),  a(k),  5(k ) 

•  final  time:  tf 
where  Vfc  6  [0 . . .  N]. 


(5.35) 

(5.36) 
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The  same  cost  function  is  used  as  the  previous  aircraft  model  with  the  addition 
of  the  controls  a(k)  and  5(k). 


N 


min 

(p(k),a(k),S(k) 


J  =  tf+p^2<p2(k) 

k= 0 


(5.38) 


The  constraints  are  also  the  same  as  Section  5.1  with  a  few  additions  including  a 
constraint  on  the  force  balance  in  the  vertical  plane,  equation  (5.45),  minimum  and 
maximum  speed,  angle  of  attack,  and  throttle  setting.  Also  the  defect  equations 
(5.39)  refer  to  the  system  dynamic  equations  (5.28-5.31) 


^k,j 

= 

0 

(5.39) 

V 

k  e  [l ...  IV],  Vj  e  [l . . .  m] 

4>min 

< 

0(^0  5;  (frmax 

(5.40) 

Rd{k ) 

< 

\/Xrel(k)  +  Vrel(k) 

(5.41) 

Vtgt(N)  -  V(N) 

= 

o 

O 

O 

(5.42) 

%rel  (A') 

= 

Rs  cos (dq^lV)  +  90°) 

(5.43) 

Vrel{N) 

= 

Rssm(^tgt(N)  +  90°) 

(5.44) 

0 

= 

Tmax^ik)  sin  a(k)  cos  4>{k)  +  L{k )  cos  4>{k)  —  mg 

(5.45) 

Vmin 

< 

v(k)  <  vmax 

(5.46) 

®-min 

< 

ot{R)  ^  oimax 

(5.47) 

$min 

< 

^{.k~}  —  ^ max 

(5.48) 

A  limit  on  dynamic  pressure  is  not  included,  because  the  maximum  speed  used  for 
the  simulations  is  a  more  conservative  constraint.  The  maximum  speed  occurs  well 
within  the  subsonic  region,  and  in  order  to  include  a  supersonic  aircraft  model,  the 
drag  coefficient  would  need  to  be  modelled  as  a  function  of  Mach  number  in  addition 
to  a. 

Rd(k )  is  still  calculated  using  equations  (5.26-5.27),  but  now  the  radar  cross  section 
is  a  function  of  xrei(k),  yrei(k),  4>(k),  T(£;),  and  a[k).  This  function  is  described  in 
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Section  2.2.4. 


5.2.3  Interceptor  Variable  Speed  and  Altitude  Model 

The  intercept  problem  is  not  limited  to  two-dimensions.  In  many  cases  the  interceptor 
should  have  at  least  a  small  altitude  offset  from  the  target,  which  can  improve  the 
view  of  certain  aircraft  features  such  as  the  horizontal  tail.  The  following  dynamics 
are  the  same  as  the  model  in  [23]  with  a  small  change  to  make  the  x  and  y  states 
with  respect  to  the  target’s  position.  Also,  the  dynamics  do  not  include  fuel  burn 
because  the  duration  of  the  trajectory  is  short.  The  same  assumptions  are  made  as 
the  previous  section  for  the  target  aircraft. 


xrei  =  v  cosy  cos  T  -  vtgt  cos  ^tgt 

ijrei  =  v  cos  7  sin  T  -  vtgt  sin  T tgt 

Tmax5  cos  a  -  D  -  yg  sin  7 

v  =  - 

d 

Tmax5  sin  a  cos  <\>  +  L  cos  4>  —  yg  cos  7 

7  =  - 

yv 

^  Tmax5  sin  a  sin  (j)  +  L  sin  (ft 
yv cos  7 

h  —  v  sin  7 


%gt  =  d 


(5.49) 

(5.50) 

(5.51) 

(5.52) 

(5.53) 

(5.54) 

(5.55) 


where  h  is  the  interceptor’s  altitude,  and  7  the  interceptor’s  flight  path  angle.  The 
same  models  for  maximum  thrust,  lift,  and  drag  are  used  as  the  previous  section. 


67 


The  variables  for  the  nonlinear  program  consist  of 


•  states  (m  =  6):  xrel(k),  yrei(k),  v(k),  7 (k),  T(/c),  h(k) 

•  controls  (n  =  3):  <f)(k),  a(k),  S(k) 

•  final  time:  tf 
where  Vfc  6  [0  ...  N], 

The  same  objective  function  is  used  for  this  model 


min  J 

(f)(k),a(k),5(k) 


N 

k= 0 


subject  to  the  following  constraints 


(5.56) 


Akj 

= 

0 

(5.57) 

V 

k  G  [1 ...  At],  Vj  G  [1 . . .  m] 

< 

—  4*ma  x 

(5.58) 

Rd(k) 

< 

\/Xrel(k)  +Vrel(k)  +  Kel(k) 

(5.59) 

Vtgt(N)  -  V(N) 

= 

O 

O 

(5.60) 

%rel  (At) 

= 

Rs  coa{^tgt{N)  +  90°) 

(5.61) 

Vrel(N) 

= 

R88m(%gt(N)  +  90°) 

(5.62) 

V min 

< 

v  ( k )  ^  vmax 

(5.63) 

&min 

< 

Oiiji)  ft  Oimax 

(5.64) 

3min 

< 

$(^)  ^  3max 

(5.65) 

'Imin 

< 

7 (k)  <  Imax 

(5.66) 

where  equations  (5.57)  are  the  defects  which  approximate  the  system  dynamics  of 
equations  (5.49-5.54). 

Now,  the  radar  detection  constraint  equation  (5.59)  includes  the  relative  altitude 
between  the  interceptor  and  target  hrei(k )  for  calculating  range.  The  range  at  which 
the  target  can  detect  the  interceptor  Rd(k )  is  also  modelled  differently 


Rd(k) 

ua 4 (k) 

1  +  e-aX'(fc) 

(5.67) 

X\k) 

=  v(k)~C(k) 

(5.68) 

v(k) 

=  xrat(k)  tan(RC A) 

(5.69) 

C  (k) 

=  \/ Keiik)  +  y2rot{k) 

(5.70) 

•E rotik ) 

=  cos(—ty tar(k))xrei(k)  -  sin tar(k))yrei(k) 

(5.71) 

l-Jrot  (  k  ) 

=  sin  (-^tar(k))xrei(k)  +  cos  tar(k))yrei(k) 

(5.72) 

As  stated  before,  the  sigmoid  form  of  equation  (5.67),  allows  the  radar  constraint 
to  be  switched  on  and  off  depending  on  whether  the  interceptor  is  inside  or  outside 
of  the  target’s  radar  cone.  cr(k)  is  a  function  of  xrei(k),  yrei(k),  hrei(k),  <t>(k),  ^(k), 
and  a(k).  A  set  of  coordinate  transformations  and  a  few  equations  in  Section  2.2.4 
map  these  variables  to  the  radar  cross  section.  The  purpose  of  equations  (5.68-5.72) 
is  to  set  X'  as  a  negative  number  if  the  interceptor  is  within  the  target’s  radar  cone, 
otherwise,  X'  is  positive.  A  graphical  representation  of  the  variables  used  to  determine 
X'  is  shown  is  Figure  5-5. 
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Figure  5-5:  (. xrot,yrot )  are  ( xrei,yrei )  rotated  by  the  negative  of  the  target’s  heading 
angle  'i’t.gt-  RCA  is  the  target’s  radar  cone  angle,  and  rj  is  the  radius  of  the  radar 
cone  at  xrot.  hrei  is  the  interceptor’s  altitude  relative  to  the  target’s.  £  is  the  distance 
from  the  target  to  the  interceptor  in  the  plane  containing  altitude  and  the  rotated 
relative  y  position  yrot. 
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5.3  Results 


5.3.1  Interceptor  Constant  Speed  Model  Results 

For  the  following  scenarios,  the  interceptor  and  target  maintain  a  constant  speed  and 
altitude.  The  interceptor’s  objective  is  to  complete  the  maneuver  in  minimum  time 
with  a  small  weight  on  bank  angle 

N 

mm  J  =  tf  +  p  4>2(k)  (5.73) 

m  tit 

Additionally,  the  interceptor  must  avoid  radar  detection  throughout  the  trajectory. 
Table  5.1  shows  the  constants  used  for  the  first  two  simulations.  In  Figure  5-6,  a 
trajectory  is  shown  where  the  interceptor  ignores  the  radar  of  the  target,  and  in 
Figure  5-7  radar  constraints  are  added  with  the  target’s  radar  detection  range  set  to 
12.2 km  when  the  interceptor’s  radar  cross  section  is  lm2.  In  exchange  for  attempting 
to  avoid  radar  detection,  the  trajectory  duration  is  11.6  seconds  longer. 


Table  5.1: 
model 


Constants  for  direct  collocation  simulations  with  constant  speed  aircraft 


N 

50 

a 

9.8  m/ s2 

0 min 

-60° 

(pmax 

60° 

*(*o) 

180° 

*tgt{to) 

0° 

V 

152.5m/s 

Vtgt 

152.5m/s 

%rel  (t0) 

27. 5km. 

Vrel  (to) 

0  km 

Rs 

3.05  km 

RCA 

45° 

P 

0.1 
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(A) 


Figure  5-6:  For  this  scenario,  the  detection  range  of  the  target’s  radar  is  set  to  zero, 
Rd(c t  =  lm2)  =  0 km,  and  the  target  is  turning  at  a  constant  rate  of  0.3°/s.  Figure 
(A)  shows  the  trajectories  of  both  aircraft,  and  Figure  (B)  is  the  bank  angle  history 
of  the  interceptor.  The  interceptor  initially  makes  a  small  turn  away  from  the  target 
followed  by  a  very  gradual  turn  towards  the  target.  At  the  end  of  the  trajectory,  a 
sharp  turn  is  applied  to  achieve  the  final  desired  relative  heading  and  position. 
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(A) 


(B) 


Figure  5-7:  In  this  case,  the  detection  range  of  the  target’s  radar  is  12.2 km  when 
the  interceptor  presents  a  radar  cross  section  of  1  m2,  and  the  target  is  turning  at 
a  constant  rate  of  0.3°/s.  Figure  (A)  shows  the  trajectories  of  the  aircraft,  and 
Figure  (B)  is  the  bank  angle  history  of  the  interceptor.  The  interceptor  begins  the 
maneuver  with  a  turn  away  from  the  target  to  fly  to  the  point,  shown  with  a  solid 
black  diamond,  at  which  it  will  cross  the  target’s  radar  cone.  Jnst  before  reaching 
this  point,  maximum  negative  bank  angle  is  applied  to  return  the  interceptor  to  a 
heading  which  will  allow  completion  of  the  maneuver.  The  interceptor  travels  along  a 
nearly  straight  path  and  then  finally  applies  maximum  positive  bank  angle  to  obtain 
the  desired  endstate. 
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Table  5.2:  Results  for  constant  speed  trajectories  including  radar  avoidance  con¬ 
straints.  The  term  Rd(cr  =  1  m2)  refers  to  the  detection  range  of  the  target’s  radar 
when  the  interceptor’s  radar  cross  section  is  1  m2. 


path 

length  (km) 

trajectory 
time  (s) 

computation 
time  (s) 

Rd(cr  =  lm2)  =  0  km 

13.6 

89.1 

15.1 

Rd(<J  =  lm2)  =  12.2  km 

15.4 

100.7 

52.6 

Figure  5-8:  Limitation  of  radar  avoidance  constraint  for  the  second  simulation, 
Rd(cr  =  lm2)  =  12.2 km.  The  interceptor  is  detected  by  the  target’s  radar  seconds 
before  crossing  the  target’s  radar  cone. 


Because  the  radar  detection  constraint  is  only  applied  at  the  nodes,  for  nearly 
every  simulation  the  interceptor  is  detected  by  the  target’s  radar  for  a  brief  moment 
before  crossing  the  radar  cone,  as  shown  in  Figure  5-8. 

There  are  a  few  variations  to  the  problem  formulation  which  might  prevent  radar 
detection  while  the  interceptor  crosses  the  target’s  radar  cone.  The  density  of  the 
nodes  around  this  area  could  be  increased,  using  the  previously  calculated  trajectory 
as  an  initial  gness  in  forming  a  new  solution.  Also,  the  optimization  problem  could 
be  solved,  increasing  the  target’s  radar  power  to  more  than  its  actual  value.  Finally, 
the  sigmoid  equation  could  be  altered  with  the  term  b  in  order  to  shift  the  radar  cone 
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Figure  5-9:  This  figure  shows  two  locally  optimal  trajectories  where  each  was  pro¬ 
duced  using  a  different  initial  guess  for  the  NLP  variables.  Trajectories  (A)  and  (B) 
are  completed  in  140.4s  and  94.2s  respectively.  The  solid  black  diamond  is  the  point 
at  which  the  interceptor  crosses  the  target’s  radar  cone. 


(k  74') 

1  _|_  e-aC(fc)  ^  1  +  e-o(c(fc)+6)  ^  J 

A  problem  with  casting  the  optimal  control  problem  as  a  nonlinear  program  is  the 
existence  of  multiple  local  minima.  Of  1, 140  simulations  run  using  MATLAB’s  Op¬ 
timization  Toolbox,  which  varied  with  initial  position  of  the  interceptor  with  respect 
to  the  target  and  target  radar  strength,  99.74%  appear  to  be  solutions  that  could  be 
globally  optimal.  For  the  three  cases  which  were  clearly  local  minima,  the  simulation 
was  performed  again,  changing  the  guess  of  the  initial  variables  at  the  nodes  by  mul¬ 
tiplying  by  1.1,  and  a  better  solution  was  obtained.  One  example  is  shown  in  Figure 
5-9. 

Another  issue  with  nonlinear  programming  is  the  time  duration  for  the  optimiza¬ 
tion  software  to  converge  to  a  solution.  Of  the  1, 140  simulations,  the  average  com¬ 
putation  time  to  reach  a  locally  optimal  solution  was  24.25  seconds  with  a  standard 
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Figure  5-10:  Impact  of  N  on  computation  time  with  the  same  constants  as  the  scenario 
in  Figure  5-7 

deviation  of  19.04  seconds.  In  compiled  implementation  and  on  a  faster  processor, 
these  computation  times  may  be  acceptable  for  real-time  implementation.  However, 
the  maximum  computation  time  was  over  6  minutes.  This  uncertainty  in  reaching  a 
locally  optimal  solution  clearly  makes  online  implementation  difficult. 

Generally,  the  computation  time  will  increase  as  the  number  of  variables  increases 
for  the  NLP.  An  important  parameter,  which  the  user  can  select,  is  N.  This  number 
represents  the  number  of  nodes  for  which  the  problem  is  discretized.  As  previously 
stated,  the  number  of  variables  for  the  NLP  is  ( N  +  1)  *  m  +  (N  +  1)  *  n  +  1  where 
m  and  n  are  the  number  of  states  and  controls  at  each  node  respectively.  Choosing 
N  is  a  balancing  process.  With  N  too  small,  large  portions  of  the  trajectory  ignore 
state  constraints  such  as  radar  detection  avoidance.  With  larger  values  of  N,  the 
computation  time  generally  increases  for  the  NLP  optimization  software  as  shown  in 
Figure  5-10.  For  the  previous  simulations  N  was  set  to  50,  in  an  attempt  to  balance 
the  computation  time  with  the  amount  of  time  that  the  interceptor  might  be  detected 
by  the  target’s  radar.  In  the  following  results  sections,  N  is  set  to  30  in  order  to  speed 
up  the  optimization  software  for  aircraft  models  with  greater  values  of  m,  n,  and  the 
number  of  constraints. 
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5.3.2  Interceptor  Variable  Speed  Model  Results 

For  the  simulations  in  this  section,  the  interceptor  is  able  to  vary  its  speed  throughout 
the  trajectory,  and  both  aircraft  are  flying  at  a  constant  altitude  of  6.1  km.  The 
interceptor  is  limited  to  flight  within  the  subsonic  region  in  order  that  the  drag 
coefficient  can  be  expressed  solely  in  terms  of  angle  of  attack.  Also,  the  interceptor 
can  only  perform  a  2-g  turn  which  is  equivalent  to  a  maximum  bank  angle  of  60°. 
With  larger  values  of  maximum  bank  angle,  the  accuracy  of  the  collocation  method 
for  approximating  the  aircraft  dynamics  decreased.  The  combination  of  speed  limits 
and  the  limit  on  bank  angle  produces  a  curve  of  maximum  turn  rates  shown  in  Figure 
5-11. 


Table  5.3:  Constants  for  variable  speed  aircraft  model  simulations  using  direct  collo¬ 
cation  _ 


N 

30 

9 

9.8m/s2 

4*min 

-60° 

4>max 

60° 

^ min 

122  m/s 

V max 

213.5m/s 

C^min 

0° 

CX-max 

30° 

$min 

0 

$max 

1 

*(*o) 

180° 

^ tgt(to ) 

0° 

v(to) 

152.5m/s 

vtgt 

152.5m/s 

%rel  {to) 

27.5  km 

Vrel  (t0) 

0  km 

Rs 

3.05  km 

RCA 

45° 

h 

6.1  km 

Ps 

0.9048/cg/m3 

P 

0.1 

mass 

9,  299 kg 

T 

±  max 

91.13kN 

The  objective  is  the  same  as  the  previous  section  where  the  interceptor  flies  the 
trajectory  in  minimum  time  with  a  small  weight  on  bank  angle 


min  J 

(f)(k),a(k),5(k) 


N 

*/+pZ>2(*) 

k=\ 


(5.75) 


Figures  5-12  and  5-14  show  two  examples  of  the  beam  intercept  and  the  corre¬ 
sponding  bank  angle  history.  For  the  first  scenario,  the  interceptor  performs  the 
maneuver  ignoring  radar  detection  constraints.  In  the  second  scenario,  the  intercep- 
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Figure  5-11:  Interceptor  maximum  turn  rate  for  variable  speed,  constant  altitude 
model 


tor  attempts  to  avoid  the  target’s  radar,  which  has  a  detection  range  of  12.2 km  when 
the  interceptor  presents  a  radar  cross  section  of  1  m2.  In  both  cases,  the  target  is 
turning  slowly  at  a  constant  rate  of  0.3°/s.  The  path  length,  trajectory  time,  and 
computation  time  for  both  scenarios  are  shown  in  Table  5.4. 


Table  5.4:  Results  for  variable  speed  trajectories  shown  in  Figures  5-12  and  5-14. 
The  term  Rd(cr  =  lm2)  refers  to  the  detection  range  of  the  target’s  radar  when  the 
interceptor’s  radar  cross  section  is  lm2. 


path 

length  (km) 

trajectory 
time  (s) 

computation 
time  (s) 

Rd(cr  =  lm2)  =  0  km 

15.6 

76.7 

175.5 

Rd(cr  =  lm2)  =  12.2  km 

16.8 

81.0 

299.5 
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Figure  5-12:  In  this  scenario,  the  interceptor  performs  the  beam  intercept  without 
considering  the  target’s  radar.  The  trajectory,  shown  in  Figure  (A),  consists  of  a 
slow  turn  away  from  the  target,  a  slow  turn  towards  the  target,  and  a  sharp  turn 
at  maximum  bank  angle  to  complete  the  maneuver.  The  interceptor’s  bank  angle 
history  is  shown  in  Figure  (B). 
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Figure  5-13:  Additional  control  and  state  histories  for  the  trajectory  in  Figure  5-12. 
The  interceptor  initially  applies  full  throttle  to  reach  its  maximum  speed.  Although 
the  maximum  speed  limit  is  met  at  the  nodes  used  in  the  NLP,  the  constraint  is 
violated  slightly  between  the  nodes.  For  the  final  30  seconds  of  the  trajectory,  the 
interceptor  slows  its  speed  in  order  to  obtain  a  faster  turn  rate. 
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Figure  5-14:  At  the  beginning  of  the  trajectory,  the  interceptor  changes  its  heading  in 
order  to  cross  the  target’s  radar  cone  at  the  point  marked  by  the  solid  diamond.  At 
this  range,  along  with  the  corresponding  radar  cross  section,  the  target  is  unable  to 
detect  the  interceptor.  Prior  to  reaching  the  black  diamond,  the  interceptor  applies 
maximum  negative  bank  angle  to  wrap  around  the  target’s  radar  cone.  Near  the  end 
of  the  trajectory,  the  interceptor  applies  maximum  bank  angle  to  obtain  the  desired 
endstate. 
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(B) 


(D) 


Figure  5-15:  Additional  control  and  state  histories  for  the  trajectory  in  Figure  5-14. 
The  interceptor  begins  the  maneuver,  applying  full  throttle  to  increase  its  speed  to 
the  maximum  limit.  Before  crossing  the  target’s  radar  cone  which  is  shown  with 
the  black  diamond,  the  interceptor  decreases  its  speed  to  obtain  a  faster  turn  rate. 
Soon  afterwards,  the  interceptor  increases  its  speed  again  to  the  limit.  During  the 
final  seconds  of  the  maneuver,  the  throttle  is  dropped  to  zero,  decreasing  the  speed 
slightly. 
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Figures  5-13  and  5-15  show  that  for  both  trajectories  the  maximum  speed  con¬ 
straint  is  violated.  This  problem  occurs  because  the  state  constraint  involving  speed  is 
only  applied  at  the  nodes.  A  more  dense  discretization  with  additional  nodes  placed 
in  the  vicinity  of  where  the  violation  occurs  could  fix  the  problem.  However,  for 
these  simulations,  the  maximum  speed  is  set  conservatively  so  that  the  aircraft  is  still 
within  the  subsonic  region  if  the  constraint  is  violated  by  a  small  amount. 

Another  potential  problem  is  the  rapid  variations  in  throttle  command.  This  could 
be  addressed  by  adding  additional  inequality  constraints  to  limit  the  rate  of  change 
between  nodes,  which  might  also  indirectly  fix  the  maximum  speed  violations. 

As  stated  in  the  previous  section,  a  shortcoming  to  creating  trajectories  with 
direct  collocation  is  the  possibility  of  converging  to  a  local  minimum.  The  maneuvers 
in  Figures  5-12  and  5-14  appear  to  be  candidates  for  the  globally  optimal  solution. 
However,  in  5.2%  of  522  conducted  simulations,  where  the  boundary  conditions  and 
radar  strength  were  varied,  the  converged  solution  is  definitely  not  a  global  minimum. 
A  better  solution  can  often  be  obtained  by  simply  varying  the  initial  guess  for  the 
NLP  as  shown  in  Figure  5-16. 

Also,  the  computation  time  to  reach  a  converged  solution  is  still  a  problem.  Of 
the  522  simulations,  the  average  computation  time  was  208  seconds  with  a  standard 
deviation  of  147  seconds.  The  maximum  time  to  reach  a  locally  optimal  solution  was 
over  41  minutes.  With  computation  times  that  are  potentially  this  long,  real-time 
implementation  seems  questionable. 
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Figure  5-16:  For  the  above  trajectories,  the  same  constants  are  used  as  in  Table 
5.4,  but  xrei(t0 )  =  28.6 km  and  yrei(to )  =  7.3 km.  The  target’s  radar  power  is  set  to 
detect  an  aircraft  at  12.2 km  when  the  radar  cross  section  is  lm2.  Figure  (A)  shows  a 
trajectory  that  obeys  all  of  the  NLP  constraints,  but  is  clearly  just  a  locally  optimal 
solution.  Figure  (B)  is  a  new  trajectory  with  the  same  initial  conditions  as  Figure 
(A),  but  the  initial  guess  for  the  NLP  variables  was  perturbed  by  110%.  This  new 
trajectory  resembles  what  could  be  the  optimal  solution. 
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5.3.3  Interceptor  Variable  Speed  and  Altitude  Model  Results 

For  this  scenario,  the  interceptor  must  perform  a  beam  intercept  in  three  dimensions. 
The  trajectory’s  initial  conditions  are  shown  in  Table  5.5.  In  addition  to  finishing 
the  trajectory  with  the  desired  relative  heading  and  position  in  the  (x,y)  plane,  the 
interceptor  must  finish  the  maneuver  with  610m  of  altitude  separation  above  the 
target.  Throughout  the  duration  of  the  maneuver,  the  target  varies  its  heading  with 
a  constant  turn  rate  and  maintains  a  constant  speed.  For  the  first  scenario  presented 

Table  5.5:  Constants  for  direct  collocation  simulations  with  variable  speed  and  alti¬ 
tude  aircraft  model _ 


N 

30 

9 

9.8  m/s2 

4*min 

1 

o 

o 

4*  max 

60° 

V min 

122  m/s 

Umax 

213.5m/s 

&min 

0° 

& max 

30° 

$min 

0 

Omax 

1 

*“) min 

-20° 

^7 max 

20° 

*(*o) 

180° 

^tgM 

0° 

v(t0) 

152.5m/s 

vtgt 

152.5m/s 

%rel  {to) 

27.5  km 

1/ relit  o) 

0  km 

R,s 

3.05  km 

RCA 

45° 

Ps 

0.9048 /cg/m3 

h(t0) 

6.1/cm 

h(tf) 

6.405/cm 

htgt 

5.795 km 

RCA 

45° 

V 

0.1 

mass 

9,  299 kg 

T 

max 

91.13/c  IV 

in  this  section,  the  target  is  turning  at  a  rate  of  0.3°/s,  and  the  interceptor  performs 
the  trajectory  without  considering  the  radar  of  the  target.  In  the  second  scenario, 
the  target  is  flying  along  a  straight  path,  and  the  interceptor  must  avoid  the  target’s 
radar  which  can  detect  an  aircraft  at  9.2 km  when  the  radar  cross  section  is  lm2. 
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Figure  5-17:  Figure  (A)  shows  a  beam  intercept  trajectory  in  the  (x,y)  plane  where 
the  target  is  turning  at  a  constant  rate  of  0.3°/s.  The  interceptor  begins  with  a  slight 
turn  away  from  the  target  and  travels  along  an  approximately  straight  line  for  the 
majority  of  the  maneuver.  The  trajectory  is  completed  with  a  final  turn  towards  the 
target.  Figure  (B)  is  the  interceptor’s  bank  angle  history. 


(A)  (B)  (C) 


Figure  5-18:  This  figure  shows  the  additional  states  and  controls  of  the  trajectory 
found  in  Figure  5-17.  Initially,  the  interceptor  climbs  to  an  altitude  significantly 
higher  than  the  desired  final  value.  It  seems  as  if  this  excess  climb  would  only  slow 
the  completion  of  the  maneuver.  This  lack  of  intuitiveness  in  the  trajectory  might 
indicate  that  the  solution  is  not  close  to  the  global  optimum. 


87 


(B) 


Figure  5-19:  Figure  (A)  shows  the  trajectory  in  the  (x,y)  plane  and  Figure  (B)  is  the 
bank  angle  history  of  the  interceptor.  The  interceptor  begins  the  maneuver  with  a 
turn  away  from  the  target  to  fly  to  a  point  along  the  radar  cone,  marked  by  a  black 
diamond,  where  it  will  not  be  detected.  Upon  reaching  this  point,  a  sharp  turn  is 
applied  to  point  the  interceptor’s  velocity  vector  more  towards  the  target.  Then,  the 
aircraft  flies  along  a  nearly  straight  path  and  finishes  the  trajectory  with  a  sharp  turn. 


Figure  5-20:  Figures  (B)  and  (E)  show  that  the  interceptor  lowers  its  speed  when  a 
large  turn  is  initiated.  The  altitude  history  shown  in  Figure  (A)  is  less  intuitive  and 
may  be  the  result  of  the  solution  having  converged  to  a  local  minimum. 


Table  5.6:  Results  for  variable  speed  and  altitude  trajectories  shown  in  Figures  5-17 
and  5-19.  The  term  Rd{(J  =  lm2)  refers  to  the  detection  range  of  the  target’s  radar 
when  the  interceptor’s  radar  cross  section  is  lm2. 


path 

length  (km) 

trajectory 
time  (s) 

computation 
time  (s) 

T tgt  =  —0.3  °/s,Rd(cr  =  lm2)  =  0  km 

15.74 

76.3 

712.9 

T  tgt  =  0.0  °/s,Rd(cr  =  lm2)  =  9.2  km 

17.55 

88.2 

659.1 

Optimal  solutions  for  this  aircraft  model  where  significantly  more  difficult  to  ob¬ 
tain  than  with  the  constant  altitude  models,  as  is  shown  by  the  computation  times 
in  Table  5.6.  There  is  some  uncertainty  as  to  whether  the  results  shown  above  are 
close  to  the  globally  optimal  solution,  because  the  altitude  history  shown  in  both  sce¬ 
narios  lacks  intuition.  For  many  simulations,  the  optimization  software  converged  to 
a  peculiar  solution  where  the  interceptor  flies  a  circuitous  route.  This  occurred  with 
larger  values  of  the  target’s  radar  strength.  In  some  cases,  providing  the  optimization 
software  with  a  better  initial  guess  seemed  to  improve  the  converged  solution.  For 
instance,  an  initial  guess  that  is  dynamically  feasible  and  meets  the  final  endstate 
constraints  usually  helped  convergence  for  the  problem  with  added  radar  constraints. 

5.3.4  Discussion  of  Results 

The  constant  speed  and  altitude  aircraft  model  simulations  had  the  fastest  conver¬ 
gence  times  and  also  had  the  highest  rate  of  converging  to  what  appears  to  be  the 
globally  optimal  solution.  For  a  real-time  implementation,  a  combination  of  the  dif¬ 
ferent  models  could  be  used.  Initially,  a  trajectory  would  be  produced  for  the  constant 
speed  aircraft  model.  Once  the  computer  has  reached  a  solution,  additional  trajecto¬ 
ries  could  be  produced  using  the  models  with  greater  degrees-of-freedom.  With  this 
approach,  the  interceptor  might  not  obtain  the  optimal  trajectory,  but  the  chances 
of  obtaining  a  feasible  trajectory  are  increased. 
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Chapter  6 


Trajectory  Interpolation 


In  this  chapter,  the  trajectory  interpolation  method  of  [12]  is  applied  for  producing 
maneuvers  which  resemble  a  beam  intercept.  The  idea  behind  this  approach  is  as 
follows.  Initially,  a  set  of  trajectories  for  a  particular  type  of  maneuver  are  produced 
offline  using  nonlinear  programming  software.  Each  trajectory  in  the  set  accomplishes 
a  similar  goal  and  closely  resembles  the  others.  The  trajectories  are  unique  in  that 
they  have  different  parameters  such  as  boundary  conditions  encompassing  a  variety 
of  values.  For  example,  a  set  of  beam  intercepts  might  be  formulated  with  each 
trajectory  containing  a  different  initial  range  between  the  target  and  the  interceptor. 
The  second  step  of  the  method  is  the  online  calculation  of  trajectories.  The  method 
interpolates  between  the  previously  calculated  trajectories  to  find  a  trajectory  for  the 
desired  parameters. 

Trajectory  interpolation  is  essentially  a  relaxation  of  a  nonlinear  parametric  pro¬ 
gramming  problem  where  the  objective  function  is  removed  [12],  Feasible  trajectories 
are  synthesized  as  an  analytic  continuation  problem  from  a  compact  parametric  rep¬ 
resentation  of  a  particular  set  of  trajectories.  In  creating  new  trajectories,  the  focus  is 
on  maintaining  feasibility.  Optimality  is  not  explicitly  considered;  however,  in  certain 
cases  the  interpolated  trajectories  are  close  to  the  optimal  solution  [13]. 
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The  set  of  feasible  trajectories  is  described  by 


h(p,  a)  =  0 

(6.1) 

g(p)  <  0 

(6.2) 

where  p  is  the  set  of  analytic  variables  which  define  a  trajectory.  In  this  thesis,  the 
values  of  p  are  b-spline  coefficients,  a  is  a  vector  which  allows  the  user  to  chose 
specific  values  for  the  parameters  which  correspond  to  a  desired  trajectory.  The 
equality  constraints  typically  describe  equations  encompassing  aircraft  dynamics  and 
boundary  conditions.  The  inequality  constraints  impose  limits  on  the  states  and 
controls  of  the  vehicle. 

Given  two  trajectories  of  similar  form,  V\  and  v2  where  v  =  [p;a],  the  method 
interpolates  from  one  to  the  other,  creating  a  similar  set  of  trajectories.  Beginning  at 
trajectory  V\ ,  the  method  interpolates  towards  v2  with  a  series  of  steps.  The  direction 
of  the  step  is  formed  by  creating  u(s)  where 

u(s)  —  v2  -  vi  (6.3) 

and  projecting  this  vector  back  onto  the  set  of  feasible  trajectories  defined  by  equa¬ 
tions  6.1  and  6.2.  The  projection  is  defined  as 

u  =  Z+u  =  Z(ZtZ)-1Zt(v2  -  m)  (6.4) 

where  Z  is  a  basis  for  the  nullspace  of  Dv [h:  g,j0] .  Dv[h;gj0]  is  a  matrix  consisting 
of  the  partial  derivatives  of  each  equality  constraint  and  active  inequality  constraint, 
which  are  denoted  with  the  set  Jo,  with  respect  to  each  of  the  variables  in  v.  Along  the 
path  towards  v2,  the  method  will  reach  the  desired  trajectory  v(i  which  contains  the 
desired  parameters  in  a.  The  method  is  shown  in  Figure  6-1  and  succinctly  described 
with  the  following  Pseudo-code  found  in  [12]. 
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Figure  6-1:  This  figure  found  in  [12]  describes  the  trajectory  interpolation  method. 
The  algorithm  begins  at  trajectory  iq  and  steps  towards  trajectory  v2-  At  each 
step  feasibility  is  maintained  by  projecting  the  direction  vector  u(s)  back  onto  the 
hyperplane  of  equality  constraints  h(p,  a)  and  active  inequality  constraints  gj0(p )  =  0. 
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Pseudo-code  for  Interpolation  Algorithm 

•  Obtain  feasible  maneuvers  v\  =  [ppan]  and  v2  =  [^2;«2]  of  the  same  type  but 
satisfying  numerically  different  boundary  conditions. 

•  Integrate  v(s)  =  u(s,v )  with  independent  variable  s  along  [0 ,agoai  —  aq]  with 
initial  condition  u(0)  =  V\  and  u(s,v )  at  any  point  s  given  by  the  following 
logic: 


1: 

u(s 

)  =  V2(s)  -  v(s) 

2: 

evaluate  Dvh(v(s )) 

3: 

Z  = 

=  basis  of  Null(Z>t,/r) 

4: 

M0  : 

=  Proj^M 

5: 

M0  « 

s—  m0  •  ( da/ds )_1 

6: 

Jl- 

=  {iVhip)  >  o},  Vi  G 

7: 

if  Ji  =  0 

8: 

U  =  Mq 

9: 

else 

10: 

J2  —  {j  £  Ji\Dvgj  •  Mo  >  0} 

11: 

if  J2  =  0 

12: 

M  =  Mo 

13: 

else 

14: 

Z'  =  basis  of  Null  ([Dvh;  Dvgj2] ) 

15: 

m  =  Proj  z'U 

16: 

m  < —  m  •  (da/ds)-1 

17: 

end 

18: 

end 

local  equality  derivatives 

tangent  space  basis 

project  difference  using  equation  6.4 

normalize  arc  length 

test  for  active  inequalities  where  M 

is  the  number  of  inequality  constraints 

none  active 

test  lst-order  inequality  behavior 
none  active  to  lst-order 

follow  lst-order  active  inequalities 
re-project  difference 
normalize  arc  length 


6.1  Aircraft  Kinematic  Model 

The  interceptor  is  modelled  as  a  constant  speed  aircraft,  and  the  target  is  assumed  to 
fly  with  a  constant  turn  rate  throughout  the  trajectory.  Radar  avoidance  constraints 
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are  not  added  to  the  problem.  The  system  dynamics  are  described  in  a  relative  frame 
using  the  same  constant  speed  aircraft  model,  found  in  Section  5.1. 


X 

=  V  cos  T  -  Vtgt  COS  T tgt 

(6.5) 

y 

=  v  sin  T  -  vtgt  sin  T tgt 

(6.6) 

T 

ta xi{4>)g 

V 

(6.7) 

*t9t 

=  d 

(6.8) 

\E 'tgt  is  found  through  simply  integrating  equation  (6.8),  and  the  states  (x,  y,  and  T) 
are  discretized  with  a  B-spline  basis  set. 


nx 


x(t) 

i= 1 

(6.9) 

v{t) 

ny 

=  J2cyyBj,ky(rJ  :  j  +  ky) 

3= 1 
n y 

=  J2ci,*BlM(T,l:l  +  k„) 

(6.10) 

T(r) 

(6.11) 

1=1 


where  r  =  1/T,  T  is  the  time  to  complete  the  trajectory,  kx  =  ky  =  k^  =  6  describe 
the  order  of  the  splines,  and  nx  =  ny  =  riq,  =  15  are  the  number  of  spline  coefficients 
for  each  of  the  variables.  The  knot  sequence  for  x,  y,  and  T  are  dehned  as  Sx  = 
Sy  =  Sq,  =  {06,  jh,  . . . ,  l6}  where  06  and  l6  represent  six  knots  at  0  and  1. 

The  b-splines  share  the  same  order  and  knot  sequence  as  those  found  in  [12],  The 
trajectory  can  now  be  described  with  the  following  set  of  variables 

p=[fe}ai.fe»}>„  {«,*}&,  ri- 

The  model  dynamics  are  sampled  over  the  set 
Se  =  {ti,  . . . ,  rn}  =  {0,  3^,  •  •  • ,  f§,  i§,  §§},  and  the  equality  constraints  are  de- 
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fined  with  the  following  equations 

1  -  v  cos('h(ri))  +  vtgt  cos(^^nT) 

,  ,  .  bx'{rn )  -  v  cos(T(rn))  +  vtgt  cos(4 ’tgtrnT) 

h\P,  a)  = 

^y'(n)  -  usin('h(ri))  +  vtgt  sin(V tgt^T) 

\  ^y'(xn)  -  usin(T(rn))  +  vtgtsm(i >tgtrnT) 

Rather  than  including  the  dynamics  for  T  in  the  equality  constraints,  a  limit  on  T  is 
included  in  the  inequality  constraints  which  are  described  later. 

The  boundary  conditions  are  defined  as 

x(0)  -  x(t0) 

2/(0)  -  y(t0) 

T(0)  -  V(t0) 

x{l)  ~  Rs  cos (i’tgtT  +  90°) 
y(  1)  -  RsSm{^tgtT  +  90°) 

T(l)  -  ^tgtT  +  90° 

where  x(t0 ),  y(t0),  and  T(f0)  are  constants  describing  the  states  at  the  initial  time, 
and  Rs  is  the  desired  final  separation  distance  between  the  interceptor  and  target. 
The  last  three  constraints  dictate  the  final  desired  relative  position  and  heading  for 
completing  the  beam  intercept. 


The  purpose  of  the  inequality  constraints  is  to  limit  the  rate  of  change  of  the  inter¬ 
ceptor’s  heading.  These  constraints  are  sampled  over  the  interval  Sineq  =  {si, . . . ,  sm}  = 


96 


{().  _ _  1 }  and  are  defined  with 


/ 


'fr'(si) 

T 


tan((/>77 \ax)Q 
v 


\ 


9  ip)  = 


^'(Sm)  _  tail {(pmax)g 
T  V 

't/'jsi)  ,  tan  {-<j>max)g 
T  '  V 


(6.14) 


V'fr'(Sm)  _  tan {<j>max)g  I 

T  V  ) 


6.2  Interpolation  Methods 

Two  methods  are  applied  for  entering  the  parameters  into  the  equality  constraints. 
In  both  cases,  a  trajectory  must  be  produced  for  the  desired  values  of  initial  relative 
y  position,  target  speed,  and  target  turn  rate.  All  other  variables  including  the  initial 
target  heading,  interceptor’s  speed,  initial  relative  x  position,  and  initial  interceptor 
heading  are  constant.  This  method  works  under  the  assumption  that  the  interceptor 
has  sufficient  time  to  fly  to  a  state  with  the  above  constants  before  starting  the  beam 
intercept. 


6.2.1  Single  Interpolation  (SI) 

The  first  method  involves  setting  the  initial  relative  y  position  as  a  and  the  target 
speed  and  target  turn  rate  as  a  function  of  a. 


a 

=  y(t  o) 

(6.15) 

Vtgt{<x) 

—  cii  T  b\0i  T  Cio; 

(6.16) 

^tgtia) 

=  0,2  T  +  C2OC 

(6.17) 
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The  coefficients  for  vtgt(a )  and  ^tgt(a)  are  found  using  the  equations 


vtgtvi  =  «i  +  havi  +  cxa2vi  (6.18) 

vtgtV2  —  a\  +  b±aV2  +  c \0?V2  (6.19) 

vtgt.Vd  =  cq  +  b\oiVd  +  c\c?Vd  (6.20) 

Vi 'tgtvi  =02  +  b-20lVl  +  C^ol1^  (6-21) 

^tgtV2  =  «2  +  b-2CXV2  +  C2«^2  (6.22) 

d 'tgtvd  —  a2  +  b20iVd  +  C20?Vd  (6.23) 


where  (aVl,vtgtvi,  ^tgtvi),  (^v2,VtgtV2,^tgtV2),  and  (aVd,  vtgtvd,  ^tgtVd)  pertain  to  the  tra¬ 
jectories  V\,V2,  and  i>d  respectively.  Now,  the  algorithm  can  create  new  trajectories 
with  varying  initial  relative  y,  target  speed,  and  target  turn  rate  by  interpolating  over 
one  variable.  V\  and  v2  are  selected  from  the  candidate  trajectory  pairs 

'Vaijjminit o);  V tgt 

min  ?  ^ tgt  max  ),  Mv  max  (to) )  Vtgt.  max  ?  tgt  min  ) 

^c(?/ram(ffl) )  V tgt 

min')  tgt  min  ),  ve(y  max  (ffi) )  ^ tgt  max  ?  ^ -tgt  max  } 

vf(y  min  (to)  i  V tgt  max  ?  ^  tgt 

max  )>  i,9(y  max  (to  )lVtgt 

min ?  t d^min^) 

Vh(y  min  (to) i  V tgt  max  1  tgt  min  ),  Vi(y  max  (to)  i  Vt.gt 

min  1  ^  tgt 

max  ) 

which  have  the  smallest  value  of  |ci|  +  |c2|  found  in  equations  6.16  and  6.17,  because 
the  method  worked  best  when  the  functions  vtgt(a)  and  tytgt(a)  were  nearly  linear. 
The  trajectories  va,i’b, . . .  ,Vi  are  produced  offline  using  NLP  optimization  software 
to  solve  for  the  set  of  p,  subject  to  the  constraints  in  equations  (6.12-6.14)  with  an 
objective  of  minimizing  the  maneuver  time. 

6.2.2  Multiple  Interpolation  (MI) 

The  second  method  also  uses  va,  ry,  ■  ■  ■  ,vt  but  interpolates  over  each  variable  in  order 
that  a  desired  trajectory  can  be  reached  by  performing  at  most  seven  interpolations. 
If  we  are  to  represent  the  sphere  of  all  engagement  trajectories  by  y(to),  vtgt  and 
T tgt,  we  can  define  a  cube  with  tuning  boundary  parameters.  Now  assume  that  the 
library  contains  an  explicit  trajectory  solution  for  each  such  extreme  condition,  which 


Figure  6-2:  Parameter  space  for  trajectory  interpolation 


is  shown  as  the  points  on  the  cube  in  Figure  6-2.  The  method  is  described  in  Figure 
6-3. 
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#1: 


#2: 


#3: 


#4: 


(to),  vtgtrnin,  ^ tgt„ 


^i(ymax(to) ,  Vtgtmin,  ^  tgtr, 


Vch  ( y  Vd  (to) ,  Vfgt  min  ">  ^  tgl  max  ) 

Vf  (i/min  (to),  V  tgtmax  ,  ^tgtmax  )  Ve  (Vmax  (to  ) ,  vtgtmax  ,  ^ tgtma. 

•- - ► - • - - - • 

Vd2(y  Vd  (to) ,  Vtgt  max  1  ^ tgt  max  ) 

Vc(ym.in(to ),  vtgtmin,  ^ tgtmin)  Vg(ymax(to) ,  Vtgtmin,  ^ tgtmin 


Vd3  (yvd  (to),  vtgtmin  ,  *tgt„ 


'Uh(ymin(to) ,  V tgtmax  ,  ^ t  tgtr, 


My  max  (to),  V tgtmax,  ^ tgtr, 


Vd4(y  Vd  (to) ,  V tgt  max  1  ^  tgtr, 


#5: 

#6: 


Vdl  (yvd  (to),  Vtgtmin  ,  ^ tgtr, 


Vd2(y  Vd  (to) ,  Vtgt  max  ?  ^ t  tgtr, 


Vd3(y  Vd  (to),  Vtgtmin  ,  tgtr, 


Vd5  (yvd  (to),  VtgtVd  ,  ^ tgtmax  ) 

Vd4(y  Vd  (to) ,  Vtgt  max  ?  ^ ttgtn 


#7: 


Vd5  (yvd  (to),  Vtgtv  ,  ^ tgtr, 


7  t'yvmax  / 

• - K - 


Vd6  (yvd  (to),  Vtgt.v, ,  ^ tgtr, 


Vd6(y  vd  (to) ,  Vtgtv 

min  / 

- - - # 


Vd(yVd  (to),  vtgtVd  ,  ^ tgtVd  , 


Figure  6-3:  The  desired  trajectory  Vd(yVd(to),vtgtv  ,  T tgt  )  is  obtained  by  performing 
three  sets  of  interpolations.  The  first  set  ($T  —  #4)  consists  of  interpolating  between 
the  eight  given  maneuvers  to  obtain  four  trajectories  which  all  have  yVd(to).  The 
second  set  (#5  —  #6)  uses  the  four  new  trajectories  to  create  two  trajectories  which 
share  yVd(to)  and  vtgtv  .  The  final  set  (#7)  interpolates  the  previous  two  trajectories 
to  obtain  i <tgtvd- 
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6.3  Results 


This  section  is  a  comparison  between  the  two  trajectory  interpolation  methods  and 
direct  collocation.  With  the  single  interpolation  (SI)  method,  the  desired  trajectory 
is  obtained  with  one  interpolation,  whereas  the  multiple  interpolation  (MI)  method 
could  require  as  many  as  seven  interpolations.  The  eight  trajectories 
(va,  Vb,  vc,  ve,  v f ,  vg,  Vh,  and  v*)  used  in  this  section  are  defined  with  the  constants 


Table  6.1:  Constants  for  trajectory  interpolation 


yminit  o) 

0  km 

ymax  (to) 

3.05  km 

152.5  m/s 

^tgtmax 

183  m/s 

d'tot 

Lyvmin 

—0.3  °/s 

\j hot 

LyLmax 

0  °/s 

x(t0) 

27.1  km 

V 

152.5  m/s 

*(to) 

IT 

^ tgt(to ) 

0 

0 max 

7r/3 

R.s 

3.1  km 

These  trajectories  vary  with  initial  relative  y,  target  speed  and  target  turn  rate 
in  order  that  trajectory  interpolation  can  be  used  to  obtain  a  maneuver  between 

[ymin(to)i  ymax(to)]i[vtgtmini  vtgtmax\i  an(^  ^ tgtmax]' 

For  the  first  simulation,  a  trajectory  must  be  found  with  yVd(to )  =  1-53 km,  vtgtv  = 
167.8 m/s,  and  ^tgtv  =  — 0.15°/s.  These  parameters  are  the  average  of  the  minimum 
and  maximum  values  for  the  given  trajectories.  For  the  SI  method,  the  trajectories 
Vciymini^iVtgt^^tgtmJ  and  ve(ymax(t0),  vtgtmax,  ^tgtmax)  are  chosen  for  Vl  and  v2 
respectively.  Figure  6-4  shows  V\  and  v2  along  with  a  few  interpolated  trajectories 
produced  in  obtaining  the  desired  trajectory  v^. 

The  trajectories  from  all  three  methods  are  shown  in  Figure  6-5.  The  solutions 
from  the  two  trajectory  interpolation  methods  are  nearly  the  same,  while  the  direct 
collocation  approach  achieved  a  solution  that  is  slightly  faster,  as  shown  in  Table  6.3. 

Figure  6-6  is  a  comparison  in  optimality  between  the  direct  collocation  method 
and  the  SI  method.  Fitting  the  trajectory  V\  with  B-splines  resulted  in  a  small  time 
increase  in  comparison  to  the  direct  collocation  method.  For  the  desired  trajectory 
v(i  the  optimal  solution  produced  by  direct  collocation  is  more  than  a  second  faster. 

The  SI  interpolation  approach  worked  well  for  various  scenarios  when  ci  and  c2 
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Figure  6-4:  The  SI  method  interpolates  between  v1(ymin(to),Vtgtmin,^tgtmirl)  and 
V2(ymax(t0),  vtgtrnax,  V tgtmax)  to  produce  the  desired  trajectory  vd(yVd(to),vtgtVd,^tgtVd)- 
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Figure  6-5:  Comparison  of  trajectory  interpolation  methods  and  direct  collocation 
for  Simulation  #1 


Table  6.2:  Simulation  #1 


path 

length  (km) 

trajectory 
time  (s) 

computation 
time  (s) 

single  interpolation  method 

12.60 

82.6 

89.7 

multiple  interpolation  method 

12.63 

82.8 

435.1 

direct  collocation 

12.41 

81.4 

16.2 
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Figure  6-6:  Optimality  difference  between  SI  method  and  direct  collocation  for  Simu¬ 
lation  #1.  The  variations  in  optimality  from  a  smooth  curve  for  the  direct  collocation 
method  are  the  result  of  the  method  having  converged  to  a  local  minimum. 

are  small  enough  so  that  the  equations  describing  vtgt(a)  and  ^tgt(a)  are  sufficiently 
linear  to  fit  within  the  area  described  in  Figure  6-7.  The  following  simulation  is  an 
example  where  this  was  not  the  case. 

For  the  second  simulation,  the  desired  trajectory  is  defined  with  the  parameters 
yVd(t o)  =  0.03 km,  vtgtVd  =  165m/s,  and  T tgtv  =  — 0.11°/s.  The  SI  method  uses  the 
trajectories  va(ymin(t0),  vtgtmin,  ^tgtmaJ  and  vb(ymax(t0),  vtgtmax,  ^tgtmin)  as  tq  and  v2 
respectively. 

As  shown  in  Figure  6-8,  the  MI  method  produces  a  solution  that  is  slightly  slower 
than  the  direct  collocation  method.  However,  the  SI  method  produces  a  trajectory 
that  does  not  appear  to  be  close  to  an  optimal  solution. 

The  SI  method  is  able  to  maintain  feasibility  in  producing  a  maneuver  with  the 
given  parameters.  However,  the  circuitous  trajectory  is  the  result  of  the  equation  for 
vtgt(a),  which  is  plotted  in  Figure  6-9.  In  order  to  reach  (ceVd,vtgtvd),  the  method  must 
interpolate  along  a  curve  which  has  a  large  value  of  |ci|.  If  the  interpolation  were 
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Figure  6-7:  The  white  region  of  the  above  contour  plot  contains  the  values  of 

[avd,vtgtvd\  for  which  vtgtmin  <  vtgt(a)  <  vtgtmax. 


Figure  6-8:  Comparison  of  trajectory  interpolation  methods  and  direct  collocation 
for  Simulation  ^2 
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Figure  6-9:  Curve  describing  problem  with  single  interpolation  approach 


Table  6.3:  Simulation  #2 


path 

length  (km) 

trajectory 
time  (s) 

computation 
time  (s) 

single  interpolation  method 

13.8 

90.7 

205.9 

multiple  interpolation  method 

13.6 

89.1 

189.7 

direct  collocation 

13.3 

87.2 

12.8 

continued  along  the  path  for  vtgt{ot)  towards  t>2,  the  target’s  speed  would  eventually 
reach  a  value  of  at  least  1  km/s.  At  this  speed,  the  interceptor  would  be  unable  to 
perform  the  desired  maneuver.  The  wobbly  trajectory  of  Figure  6-8  is  a  precursor  of 
interpolating  towards  infeasible  trajectories. 


6.3.1  Discussion  of  Results 

Fortunately,  in  certain  regions  where  one  method  struggles,  the  other  excels.  For 
the  second  simulation,  the  MI  method  obtained  a  solution  faster  than  the  SI  method 
despite  having  to  perform  seven  interpolations.  The  MI  approach  is  fastest  when  the 
desired  trajectory  has  parameters  which  are  close  to  the  given  trajectories’. 

As  shown  in  the  first  simulation,  the  MI  method  consumes  significantly  more  time, 
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over  400  seconds,  when  the  desired  trajectory  has  parameters  which  are  furthest  away 
from  the  given  trajectories’.  This  is  exactly  where  the  SI  method  works  well,  because 
the  equations  for  vtgt(oi)  and  \E 'tgt(ot)  are  more  linear. 

In  every  simulation  performed,  the  computation  time  for  the  direct  collocation 
method  was  significantly  faster  than  the  interpolation  methods’.  Both  the  SI  and  MI 
methods  could  be  improved  by  increasing  the  number  of  available  trajectories  to  be 
used  for  interpolation.  A  denser  set  of  trajectories  would  increase  the  region  where 
vtgtmin  <  vtgt(a)  <  vtgtrnax  and  T tgtrniri  <  ^tgM)  <  ^ tgtmax •  Also,  more  trajecto¬ 
ries  would  decrease  the  computation  time  for  the  MI  method,  because  the  desired 
parameters  might  occur  closer  to  a  given  trajectory’s. 

Another  issue  worth  considering  is  the  optimality  gap  between  the  direct  col¬ 
location  and  trajectory  interpolation  methods.  Because  there  is  no  guarantee  that 
the  converged  solutions  for  va , Vb,  ■  ■  ■  ,vt  are  global  optima,  better  solutions  for  these 
trajectories  might  reduce  the  maneuver  times  for  the  interpolated  trajectories. 
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Chapter  7 


Conclusions 


This  research  presents  three  different  approaches  for  creating  trajectories  which  re¬ 
semble  the  fighter  pilot  maneuver  frequently  referred  to  as  a  beam  intercept.  The 
tactic  is  typically  used  to  visually  identify  a  suspected  aircraft,  and  an  ideal  trajec¬ 
tory  is  completed  in  minimum  time  while  avoiding  detection  from  the  other  aircraft. 

This  chapter  includes  a  summary  of  the  thesis  and  ideas  for  future  work. 

7.1  Summary 

In  Chapter  2,  a  description  of  the  five  phases  of  air-to-air  combat  is  presented.  Tech¬ 
nology  and  tactics  have  changed  throughout  history,  but  these  phases  are  still  relevant. 
However,  the  importance  of  each  phase  possibly  has  shifted.  The  maneuver  phase  or 
dogfight  in  recent  conflicts  did  not  occur  nearly  as  often  as  in  WWI  and  WWII.  The 
detection  and  closing  phase  remain  extremely  important.  An  early  detection  provides 
the  fighter  pilot  with  an  opportunity  to  make  the  first  decisions.  Despite  radar  tech¬ 
nology  that  can  detect  certain  aircraft  hundreds  of  miles  away  and  missiles  which  can 
reach  far  beyond  a  pilot’s  visual  range,  the  closing  phase  is  still  essential.  A  fighter  pi¬ 
lot  might  be  forced  to  merge  with  an  electronically  confirmed,  enemy  aircraft  because 
of  the  rules  of  engagement. 

Presented  with  this  scenario,  the  pilot  must  be  able  to  perform  a  maneuver  to 
visually  identify  the  aircraft.  Four  different  methods  are  discussed,  and  the  beam 
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intercept  is  usually  the  maneuver  of  choice  clue  to  its  quick  execution  and  stealthy 
approach. 

In  Chapter  3,  five  different  trajectory  generation  methods  are  considered  for  cre¬ 
ating  maneuvers  which  resemble  a  beam  intercept.  Ideally,  an  analytic  solution  could 
be  produced.  This  would  allow  the  problem  to  be  solved  in  real-time;  however,  it 
is  unlikely  that  such  a  solution  exists  given  the  complexity  of  the  radar  avoidance 
constraints.  A  shooting  method  approach  is  also  considered  but  not  applied  in  this 
thesis  due  to  the  difficulty  of  including  path  constraints.  The  three  applied  methods 
in  this  thesis  are  mixed  integer  linear  programming  (MILP),  direct  collocation,  and 
trajectory  interpolation. 

In  Chapters  4-6,  the  three  chosen  methods  are  applied.  The  MILP  approach  has 
the  advantage  that  a  solution  can  be  produced  in  real-time.  However,  the  complexity 
of  the  radar  avoidance  constraint  is  limited.  Additionally,  the  aircraft  model  must  be 
cast  in  linear  form,  which  might  limit  the  accuracy  of  an  aircraft  model  with  more 
degrees  of  freedom  such  as  being  able  to  change  altitude. 

The  direct  collocation  approach  is  able  to  handle  the  limitations  that  MILP  faces. 
A  more  complex  radar  avoidance  constraint  and  aircraft  model  are  included  in  this 
formulation.  However,  because  the  problem  is  solved  as  a  nonlinear  program  the 
major  advantages  with  MILP,  such  as  nearly  guaranteed  real-time  convergence,  are 
lost.  In  many  cases,  the  optimization  software  is  able  to  converge  quickly  to  what 
appears  to  be  a  global  optimum.  However,  it  is  the  few,  faulty  cases  which  make  an 
onboard  implementation  difficult. 

The  trajectory  interpolation  approach  is  the  third  method  applied.  With  this 
formulation,  the  goal  is  not  to  find  an  optimal  trajectory.  Rather,  the  method  inter¬ 
polates  between  given  trajectories  with  varying  boundary  conditions  and  dynamics, 
obtaining  new,  feasible  maneuvers.  The  time  to  reach  a  solution  can  be  estimated, 
and  the  ability  of  the  method  to  reach  feasible  maneuvers  can  be  tested  offline. 
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7.2  Future  Work 


Direct  collocation 

For  the  direct  collocation  simulations,  the  initial  guess  for  the  nonlinear  program  can 
affect  greatly  the  converged  solution.  This  was  especially  true  for  the  variable  speed 
and  altitude  aircraft  model.  In  these  simulations,  the  method  often  struggled  to  reach 
a  feasible  solution.  Even  when  the  method  converged,  it  is  questionable  whether  or 
not  the  produced  trajectory  is  a  global  optimum. 

A  potential  solution  to  this  problem  is  to  use  a  genetic  algorithm  approach  to 
produce  initial  guesses.  Used  alone,  genetic  algorithms  are  not  competitive  with 
gradient  methods  for  producing  an  optimal  solution  [3].  However,  in  [37]  the  authors 
found  that  combining  genetic  algorithms  with  direct  collocation  can  improve  the 
chances  of  converging  to  the  global  optimum.  The  idea  is  that  the  genetic  algorithm, 
while  not  typically  able  to  End  a  locally  optimal  solution,  is  able  to  search  globally 
the  set  of  feasible  solutions.  The  best  solution  of  the  genetic  algorithm  is  then  applied 
as  the  initial  guess  for  the  nonlinear  program.  This  combined  approach  may  yield 
better  solutions  for  beam  intercept  trajectories. 


Trajectory  interpolation 

The  trajectory  interpolation  approach  could  also  be  explored  further  by  adding  radar 
avoidance  constraints  and  a  more  complicated  aircraft  model.  With  additional  con¬ 
straints  and  variables,  it  would  be  interesting  to  note  if  the  method  can  still  obtain 
feasible  trajectories. 

Another  aspect  of  trajectory  interpolation  to  be  explored  is  the  computation  time 
in  forming  new  trajectories.  In  the  presented  results,  the  method  was  usually  signifi¬ 
cantly  slower  than  direct  collocation.  Increasing  the  number  of  trajectories  calculated 
offline,  in  order  that  the  method  does  not  have  to  interpolate  as  far,  is  one  potential 
solution.  Defining  how  many  trajectories  is  appropriate  for  a  given  scenario  could  be 
an  interesting  topic  of  research. 
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Intelligent  target 


For  all  of  the  approaches  presented,  the  trajectory  of  the  target  aircraft  was  limited 
to  trim  flight  in  steady  level  or  a  constant  turn  rate.  Future  research  could  include 
a  more  sophisticated  target  which  would  take  precautions  to  avoid  a  beam  intercept. 
For  example,  the  target  could  respond  defensively  or  offensively  when  the  interceptor 
violates  the  radar  range  constraint.  It  could  also  intelligently  vary  its  heading  so  that 
its  radar  cone  does  not  produce  constant  blind  spots. 

Multiple  Aircraft 

Radar  blind  spots  can  also  be  avoided  by  flying  in  formation.  It  is  unlikely  that  the 
target  aircraft  will  always  fly  alone.  An  ideal  trajectory  would  include  constraints 
considering  the  radar  of  all  present  aircraft.  Additionally,  the  method  should  include 
which  aircraft  to  pursue  first. 

The  interceptor  will  most  likely  fly  with  support  aircraft  as  well.  The  problem 
could  include  additional  friendly  aircraft  which  would  be  able  to  assist  in  performing 
visual  identifications.  The  MILP  approach  can  easily  handle  more  aircraft.  However, 
this  might  be  difficult  to  include  with  the  direct  collocation  and  trajectory  interpola¬ 
tion  approaches. 

Additional  constraints 

Another  potential  problem  is  no  fly  zone  areas.  For  instance,  an  optimal  trajectory 
should  avoid  areas  densely  populated  with  surface-to-air  missile  sites.  If  these  con¬ 
straints  are  included,  the  relative  frame  formulation  used  for  the  direct  collocation 
and  trajectory  interpolation  approaches  may  be  inappropriate.  An  inertial  reference 
frame  would  facilitate  the  formulation  of  these  constraints. 
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Appendix  A 


Collocation  Example 


The  following  equations  are  an  example  of  the  collocation  method  described  in  section 
5.1.1.  The  system  dynamics  are  the  same  as  section  5.2.1 


Xrel 

=  v  cos  T  -  vtgt  cos  T tgt 

(A.l) 

Vrel 

=  v  sin  T  -  vtgt  sin  T tgt 

(A.2) 

T 

tan(0)gf 

V 

(A.3) 

^ -tgt. 

=  d 

(A.4) 

The  problem  is  discretized  with  N  +  1  nodes  so  that  the  NLP  is  formulated  with 
the  variables  xrei(k),  yrei(k),  \1 /(&),  and  tf  where  \/k  G  [0 ...  N], 

The  defects  are  calculated  using  the  equations  from  section  5.1.1 

xc  =  (xi  +  x2)/2  +  T(/i -/2)/8  (A. 5) 

x'c  =  -3(xi  -  x2)/2T  -  (/i  +  /2)/4  (A. 6) 

A  =  fc-x'c  (A. 7) 


where  T  =  %. 

The  first  step  in  forming  the  defects  for  this  problem  is  to  calculate  the  interpolated 
values  of  the  states  and  controls  which  are  expressed  in  the  dynamics.  In  this  case, 
they  are  T  and  (j) 
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'&(k)  +  '&(k  +  l)  ftan(j)(k)g  tan<p(k  +  l)g\ 

- 2 - +  T{— - v - )/8  (A‘8) 

<j>(k)  +  <j>(k  +  1) 
he  [o ... iv - 1] 

The  next  step  is  to  calculate  the  actual  defects  for  each  of  the  dynamics  equations. 
Collocation  is  not  required  for  the  target’s  heading  because  equation  (A. 4)  can  be 
integrated  easily. 


fc,xrel(k)  =  vcos<l>c(k)  -  VigtCOS^igt  (A.10) 

Xc,xrel(k)  =  -3 (xrel(k)  -  xrel(k  +  1))/2T  + 

-[(u  cos  T(/c)  -  vtgt  cos  ^tgt)  + 

(v cos  \l/ (k  +  1)  -  Vtgt  cosdh9t)]/4  (A.ll) 

Afc,i  =  fc,xrel(k)  -x'c^rel{k)  (A.  12) 

V  ke[0...N-l] 


*c{k)  = 

4>c(k)  = 
V 


fc,yrel(k)  =  v  sin  4/c(fc)  —  vtgt  sin  ^tgt  (A.13) 

x'c,yrel(k)  =  -3 (yrei(k)  -  yrei(k  +  1))/2T  + 

-  [(u  sin  T  (k)  -  vtgt  sin  Vtgt)  + 

(v  sin  T(A;  +  1)  -  vtgtsm^tgt)\/A  (A.  14) 

Afc,2  =  fc,yrel{k)  -x'Cyrel{k)  (A.  15) 

V  ke[0...N-l] 
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tan  4>c(k)g 


c,\I/(/c) 


<*W  =  -3(«W  **(*  +  1))/2T  + 


tan  c i>(k)g 


v 

■  x' 


+ 


tan  0(/c  +  l)g 


/4 


Afc,s  =  fc,*(k)-x'c^k) 

V  fc  e  [0 . . .  N  -  1] 


(A. 16) 

(A.  17) 
(A.18) 


Now,  the  dynamics  of  the  system  are  met  approximately  when  equations  A.  12, 
A.  15,  and  A.18  are  equal  to  zero. 
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